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13.  ABSTRACT  (Maximum 200 word*;  The  project  is  aimed  at  seismic  surveillance  as  part  of  on¬ 
going  efforts  for  improving  nuclear  test  ban  verification  capabilities.  The  problem 
is  complex  in  the  sense  that  underground  explosions  are  most  efficiently  monitored 
by  seismic  means,  but  that  the  distinction  between  signals  emitted  by  natural  earth¬ 
quakes  and  explosions  remains  unclear,  at  least  at  local  and  regional  distances.  In 
other  words,  seismic  wave  propagation  in  a  heterogeneous  Earth  may  easily  mask 
specific  source  signatures. 

In  Section  2,  we  report  on  major  results  from  a  mar:’ ne  seismic  reflection  survey, 
in  the  Skagerrak  Sea  and  the  outer  Oslo  Fjord.  The  profiling  area  "covers"  the 
Skagerrak  Graben,  the  southward  extension  of  the  Oslo  Rift  in  which  vicinity  both 
the  NORSAR  and  NORESS  arrays  are  located.  An5rway,  extensional  simple  shear  deform¬ 
ation  of  the  entire  lithosphere  has  been  postulated  by  Wernicke  (1985)  and  others, 
but  up  to  now  unequivocal  seismic  evidence  in  support  of  this  hypotheses  has  been 
lacking  (Kusznlr  and  Egan,  1989).  Here  we  give  evidence  on  well  defined  seismic 
reflectors  in  the  recordings  coda  correlations  seldom  exceed  0.4  units.  Since 
various  kinds  of  homogeneous  crust/lithosphere  models  seemingly  are  unable  to  account  for 
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prominent,  observational  coda  features  the  crust /lithosphere  system  were 
described  in  terms  of  self-similar  von  Karman  inhomogeneous  media  with 
correlation  distance  of  10  km  and  velocity  perturbations  of  2-4  per  cent. 
Here,  Moho  undulations  were  represented  as  a  ID  von  Karman  function. 
Synthetics  computed  for  such  kinds  of  inhomogeneous  media  quantitatively 
accounts  for  major  coda  features  like  relative  high  amplitudes,  long 
duration  and  poor  coherency  beyond  a  sensor  spacing  of  a  few  kilometres . 
Also,  P-to-S  or  S-to-P  conversions  in  the  source  vicinity  are  relative 
efficient  for  inhomogeneous  media  as  compared  to  homogenous  media.  The 
other  work  presented  in  Section  4  deals  with  fundamental  mode  Rayleigh  (Rg) 
wave  propagation  in  the  uppermost  part  of  the  crust.  The  observational 
data  comprise  recordings  from  7  arrays  located  on  4  continents.  The 
inversion  of  the  individual  phase  velocity  observations  were  considered  for 
one  layer  and  two  layers  models  above  half space.  In  the  latter  case,  the 
layer  thicknesses  were  fixed  at  0.5  km  and  1.0  km.  Unknowns  were  layer  and 
half  space  shear  velocities  and  in  case  of  Model  1  also  layer  thickness. 

The  average  of  estimated  half space  velocities  was  3.56  kms'^  and  array 
differences  were  small  here.  Also  the  Model  1  layer  shear  velocity 
estimated  were  rather  consistent  between  the  various  arrays  and  the  average 
was  2.87  kms'^.  In  contrast  the  Model  1  layer  shear  thicknesses  varied 
considerably  with  extreme  values  of  0.12  km  for  Yellowknife  and  1.6  km  for 
Alice  Springs.  The  two  layer  model  gave  a  slightly  poorer  fit  to  the 
observations  in  comparison  to  the  one  layer  model.  In  order  to  provide  a 
better  insight  in  Rg-propagation  in  an  uppermost  crust  low  velocity  layer 
(LVL)  of  thickness  1.4  km  we  generated  2D  FD  synthetics  both  for  homogenous 
and  inhomogeneous  media.  In  the  former  case,  rather  classical  surface  wave 
trains  were  produced  including  a  prominent  Airy  phase  at  c .  1.1  sec . 
period.  For  the  inhomogeneous  media  multiple  scattering  became  pronounced 
with  S-scattering  wavelets  interfering  with  the  Rg-waves.  Also,  the  Airy 
phase  was  considerably  less  distinct  while  Rg  excitation  was  dr2unatively 
reduced  for  focal  depths  below  2-3  km. 

Our  first  attempt  on  seismic  discrimination  between  earthquakes  and 
underground  nuclear  explosions  is  documented  in  Section  5.  This  problem  is 
formulated  as  an  exercise  in  pattern  recognition  approach  analysis.  An 
advantage  of  our  procedure  is  flexibility,  by  combining  both  adaptive  noise 
suppression  and  event  classification  incorporating  feature  selection 
criteria.  The  procedure  has  been  applied  to  a  learning  set  of  44  nuclear 
explosions  (8  test  sites)  and  35  earthquakes  in  Eurasia  recorded  at  the 
NORESS  array.  The  signal  features  considered  were  the  normalized  power  in 
8  spectral  bands  in  the  0. 2-5.0  Hz  range  of  the  P-wave  (6  sec)  and  the  P- 
coda  (30  sec).  Physically,  it  means  that  we  exploit  potential  differences 
in  the  shape  of  earthquake  and  explosion  spectra,  respectively.  Other 
features  included  are  peak  P  and  P-coda  amplitude  frequencies  and  relative 
P/P-coda  power.  These  19  features  were  extracted  either  from  conventional 
array  beam  traces  or  optimum  group  filtered  traces  (OGF-removal  of  coherent 
low-frequency  noise) .  Using  the  feature  selection  algorithm,  based  on 
estimates  of  the  expected  probability  of  misclassification  (EPMC),  only  2 
to  4  features  were  needed  for  optimum  discrimination  performance.  The 
dominant  features  were  coda  excitation  and  P-  and  P-coda  power  at  lower 
signal  frequencies.  Furthermore,  feature  parameters  extracted  from  the  OGF 
traces  had  a  slightly  better  performance  in  comparison  to  those  extracted 
from  beam  traces.  Finally,  there  were  no  misclassifications  for  OGF- 
derived  features  when  the  explosions  population  was  limited  to  E.  Kazakh 
events,  while  including  events  from  the  other  test  sited  lead  to  a  decrease 
in  discrimination  power. 
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Skagerrak  Sea  which  ate  interpreted  as  major  low  angle  shear  zones  (The 
Bamble  Fault  -  BF)  also  underlying  the  Skagerrak  graben  segment  of  the  Oslo 
Rift.  A  few  of  them  can  be  traced  almost  continuously  from  mid-crust 
tnrough  the  lower  crust,  off-setting  Moho  and  continuing  downwards  to  50-60 
km  depth  (16  sTWT)  .  This  an-i  similar  results  from  1730  km  of  deep  seismic 
lines  in  the  Skagerrak  Sea  indicated  that  the  crust  and  mantle  did  inherit 
a  pronounced  structural  fabric  from  the  earlier  Proterozoic  Sveconorwegian 
orogeny.  Reactivation  of  such  weakness  zones  surfaced  during  the  rifting 
may  explain  the  occurrence  of  simple  shear  in  the  lower  crust  below 
Skagerrak . 

In  another  study,  based  on  the  same  R/V  Mobil  Search  profiling  records,  we 
focused  on  mapping  faults  in  the  upper,  presumed  brittle  part  of  the  crust. 
The  northeastern  area  of  the  Skagerrak  Sea  appeared  most  suitable  for  this 
undertaking  because  the  sediment  coverage  here  is  very  thin  or  absent. 

Close  inspections  of  our  high  quality  sections,  reprosessed  to  full 
commercial  exploration  standard,  gave  strong  evidence  of  seismic  images  of 
an  abundance  of  moderately  steep  (20  -  40  deg)  faults  in  the  crystalline 
basement.  The  fault  geometry  could  be  inferred  at  profile  intersections, 
and  moreover  was  found  to  coincide  with  the  general  tectonic  imprints 
observed  on  land  and  along  the  rim  of  the  Skagerrak  Graben.  Fault 
interspacing  was  ca.  3-5  km  and  depth  penetration  10-15  km. 

In  a  seismic  surveillance  context,  the  above  results  are  considered 
important  on  two  accounts;  firstly  an  extensively  cracked  crust  is  likely 
to  generate  considerably  amount  of  shear  wave  energy  even  for  simple  P-type 
explosive  sources.  Observational  evidence  stem  from  quarry  blast 
recordings  where  often  the  largest  signal  amplitude  is  found  on  the 
transverse  component.  This  is  attempted  synthesized  by  putting  P-sources 
in  inhomogeneous  media  (Section  4).  Secondly,  earthquake  occurrence  is 
often  related  to  faults  and  lineaments  derived  from  satellite  images  of  the 
Earth  surface.  This  could  be  a  dubious  undertaking  in  view  of  an  epicenter 
determination  accuracy  seldom  better  than  5-10  km  and  that  projection  of 
"surface"  faults  deep  into  the  crust  is  not  well  constrained  as 
demonstrated  for  the  Skagerrak  Sea. 

In  previous  work  under  this  contract  we  have  presented  novel  schemes  for 
epicenter  locations  and  signal  detections  using  3-component  observations. 

In  section  3  an  experiment  is  described  by  which  we  use  the  detector 
described  in  Ruud  and  Husebye  (1992)  for  automatically  picking  P-  and  S- 
arrivals  in  local  event  records  stemming  from  the  Norwegian  Seijmic 
Network.  For  automatic  epicenter  determination  a  robust  grid-search  method 
well  suited  for  estimation  problems  with  non-Gaussian  observational  errors 
is  introduced  in  order  to  handle  outliers.  Even  several  I'rge  arrival  time 
outliers  did  not  prevent  solutions  close  to  the  "true"  epicenter.  In  our 
experiment,  38  local  events  from  tite  August  1991  bulletin  were  located. 

The  number  of  detecting  stations  varied  from  3  to  10  out  of  a  total  of  15 
stations.  P-  and  S-picking  errors  were  small,  mostly  within  0.5  sec  for 
both  P  and  S.  Phase  identification  eriors,  causing  severely  wrong  P-  and  S- 
arrivals  were  more  frequent.  Decent  epicenter  determinations  were  obtained 
even  for  events  with  30-50  per  cent  outliers.  The  RMS  location  difference 
of  our  "automatic"  solutions  compared  to  those  i.n  the  manually  produced 
bulletin  were  about  +/-  15  km. 
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In  Section  4  we  present  two  works  addressing  the  problem  of  seismic  wave 
propagation  in  inhomogeneous  media.  Firstly,  details  are  given  on  the 
methodology  used  for  computing  synthetic  seismograms  from  2D  finite 
difference  (FD)  solutions  of  the  elastic  wave  equation.  The  principal  aim 
of  this  study  was  to  explore  the  range  of  crust/lithosphere  structural 
models  which  were  capable  of  producing  strong  and  persistent  coda  waves 
commonly  observed  at  local  and  regional  distances.  The  starting  model  was 
homogenous  but  included  a  Moho  bump  like  that  associated  with  the  Oslo  Rift 
in  the  vicinity  of  the  NORESS  siting  area.  The  corresponding  synthetics 
were  similar  to  those  obtainable  using  ray  tracing  methods  or  in  other 
words  no  significant  coda  waves  were  generated.  Introducing  crustal 
velocity  gradients  in  the  crust  neither  contributed  much  to  the  synthetic 
record  complexity.  In  contrast,  a  model  with  a  sinusoidal  shaped  Moho  (X  = 
8  km.  Amp  =  1  km)  produced  an  abundance  of  multiple  reflected  and  scatlered 
P-  and  S-phases  giving  the  appearance  of  strong  and  persistent  P-  and  S- 
wave  codas.  The  corresponding  coda  correlations  over  a  line  array  length 
of  5-10  km  were  rather  high  (around  0.5  -  0.7  units)  while  for  real 
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Report  summary 


Task  objectives; 

Very  detailed  crust/lithosphcre  mapping  by  seismic  profiling  means 

Automatic  analysis  of  seismograph  network  data  and  bulletin  production 

Seismic  wave  propagation  in  the  lithosphere  using  2D  finite  difference  solutions  of 
the  elastic  wave  equation 

Enhanced  seismic  source  discrimination  at  teleseismic  distances  using  NORESS 
array  data 


Technical  and  Scientific  Problems 

In  1987  one  of  us  (E.S.  Husebye)  had  at  his  disposal  the  seismic  ship  R/V  Mobil  Search 
(Courtesy  of  Mobil  Exploration  Inc  (Norway))  for  30  days  of  profiling  in  the  Skagerrak 
Sea.  The  profiling  area  includes  the  southern  extension  of  the  Oslo  Rift  in  which  vicinity 
the  NORESS  array  is  located.  Use  of  the  seismic  reflection  data  (totalling  1730  km)  have 
enabled  us  to  image  crust  and  lithosphere  structures  in  exceptional  detail,  including  Moho 
disruptions.  Likewise,  a  large  number  of  minor  and  major  faults  have  been  mapped  both 
within  the  crystalline  crust  and  also  in  the  sub-Moho  part  of  the  lithosphere.  Two  of  the 
report  works  deal  with  this  kind  of  seismic  mapping. 

In  previous  reports  we  have  addressed  problems  related  to  more  efficient  signal  detection 
and  also  robust  techniques  for  epicenter  locations.  This  kind  of  methods  are  essential  for 
automating  seismograph  network  operations  including  that  of  bulletin  production  without 
analyst  interference.  In  this  report  wc  demonstrated  that  the  latter  problem  is  solvable  in 
practise  using  observational  data  from  the  Norwegian  Seismograph  Network. 
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The  nature  of  seismic  wave  propagation  at  local  and  regional  distances  are  important  for 
the  optimum  design  of  seismic  source  classification  criteria.  In  this  context,  we  have 
explored  the  usefulness  of  2D  finite  difference  .solutions  of  the  elastic  wave  equation  in 
combination  with  varying  degrees  of  structural  complexities  in  the  crust  and  lithosphere. 
The  major  problem  is  to  generate  synthetics  which  at  least  approximately  match  typical 
observational  features  like  strong  coda  waves  of  extended  duration.  In  this  re.spect,  self¬ 
similar  scattering  media  of  von  Karman  type  proved  far  more  efficient  for  coda  generation 
that  classical  stratified  .media.  In  a  Rg-sludy,  ba.sed  on  recordings  from  7  arrays  on  4 
continents,  we  also  used  synthetics  for  gaining  a  better  insight  in  Rg-propagation  in 
complex  media. 

The  major  challenge  in  seismic  surveillance  in  the  amtext  of  a  potential  test  ban  treaty,  is 
that  of  classifying  events  as  cither  an  earthquake  or  man-made  explosion.  This  kind  of 
problems  were  formulated  as  an  exercise  in  pattern  recognition  approach  analysis.  Using 
NORESS  array  recordings  of  telcscismic  Eurasian  events,  our  analysis  procedure  combined 
both  adaptive  noise  suppression  and  event  cla.ssification  incorporating  feature  selection 
criteria.  The  signal  features  considered  were  the  normalized  power  in  8  sjsectral  bands  in 
the  0.2  -  5.0  Hz  range  of  the  P-wavc  (6  .sec  duration)  and  the  subsequent  P-coda  (30  sec 
duration).  Additional  features  were  peak  P  and  P-coda  amplitude  frequencies  and  relative 
P/P-coda  power.  Out  of  these  19  features  only  a  subset  of  2-4  features  were  needed  for 
optimum  di.scrimination  performance  when  tested  on  44  presumed  nuclear  explosions  and 
35  presumed  earthquakes. 

Methodology  used  in  the  seismic  data  analysis 

The  .seismic  profiling  data  u.sed  in  crustal  .structure  imaging  were  processed  using 
commercial  "oil  prospecting"  software.  The  final,  refined  data  procc.ssing  were  performed 
at  BIRPS,  Cambridge  Univ.,  UK  which  cxpcrti.se  in  the  particular  field  of  signal 
proce.s.sing  is  recognized  world-wide.  Regarding  automatic  network  operation,  a  novel, 
robust  "grid  search"  epicenter  location  scheme  has  been  introduced.  It  has  been  tested  on 
kx:al  network  data,  and  provided  decent  location  c.stimatcs  even  in  cases  of  up  to  40  per 
cent  gro.ss  errors  in  the  input  parameters  like  P-  and  S-arrival  mes  and/or  azimuth. 


Regarding  our  synthetics  seismogram  analysis,  2D  finite  difference  techniques  were  used 
by  geophysicists  more  than  two  decades  ago.  However,  a  unique  feature  of  our  approach 
is  that  the  model  sizes  used  enable  us  to  simulate  real  recording  at  local  distances  (out  to 
3-400  km).  The  event  classification  techniques  used,  albeit  reflecting  a  pattern  recognition 
approach,  proved  to  be  both  flexible  and  powerful  due  to  the  incorporation  of  feature 
selection  criteria. 
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Major  Scientific  Results 


In  Section  2  works  bearing  on  imaging  (1-3  km  resolution)  of  crust/lithosphcre  structures 
are  presented.  We  were  here  able  to  map  many  faults  in  the  upper,  brittle  crust  with  dip 
angles  in  the  range  20-40  deg.,  depth  extent  10-15  km  and  (horizontal)  spacing  of  3-5  km. 
The  most  spectacular  finding  was  the  socallcd  Bamble  Fault  which  could  be  traced 
through  the  brittle  crust,  continuing  into  the  presumed  lower  ductile  crust,  off-.sctting 
Moho  and  continuing  further  down  to  50-60  km  depth  (16  sTWl').  Note,  seemingly  only 
seismic  reflection  profiling  can  provide  such  detailed  structural  image,  definitely  not  using 
other  geophysical  means  nor  from  geological  mapping.  Clearly,  an  explosive  source  in  a 
"cracked"  upper  crust  would  generate  a  considerable  amount  of  shear  waves  thas 
quantitatively  explaining  the  observational  fact  that  for  this  type  of  sources  the  strongest 
signals  are  often  found  on  the  transverse  component. 

A  long  standing  goal  in  the  .seismological  community  is  that  of  automating  seismograph 
network  operation.  This  problem,  dealt  with  in  Section  3,  has  become  more  acute  In 
recent  years  simply  because  modern,  highly  .sensitive  seismometers  record  daily  a  very 
large  number  of  man-made  local  explosions.  Our  novel  v/eighted  "grid-search"  epicenter 
location  schemes  was  tested  on  Norwegian  Seismograph  Network  (NSN)  recording  for 
Aug.  1991.  The  input  parameters  here,  F’-  and  S-arrival  times  and  occasionally  P-azimuth 
were  extracted  from  records  using  the  Ruud  and  Hu.scbyc  (BSSA  1992)  detector.  The 
location  scheme  proved  robust  enough  to  handle  partly  faulty  input  parameters.  The  most 
common  detector  errors  were  a.ssociated  with  faulty  phase  identifications;  for  example  by 
mi.ssing  a  weak,  first  arriving  P-pha.se  the  subsequent  and  often  strong  S-phasc  would  be 
denoted  P.  Pha.se  picking  as  such  appeared  to  be  accurate:  -r/-  0.5  see  and  -h/-  1.0  sec  for 
P-  and  S-phases  respectively.  A  detail  comparison  with  analy.st  epicenter  solutions  (HYPO 
71)  for  38  events,  gave  that  the  average  difference  in  epicenter  coordinates  were  2  km  (lat, 
long)  with  a  std.  deviation  of  15  km.  Important,  with  this  kind  of  epicenter  accuracy  is 
would  be  easy  to  introduce  advanced  signal  prtxxssing  schemes  in  the  post-detection  state 
simply  bccau.se  potential  pha.se  arrivals  like  Pn,  Pg.  PmP,  pP,  Sn,  Sg,  SmS  etc  can  be 
accurately  predicted  within  lime  frames  of  at  most  a  few  seconds. 
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In  Section  4  two  work  bearing  on  synthesizing  crust/litliosphere  seismic  wave  propagation 
are  described.  A  basic  problem  is  strong  P-  and  S-codas,  a  typical  observational  feature  at 
local  distances,  arc  not  reproducible  using  conventional  stratified  cru.st/lithosphere  models. 
For  our  synthetics  calculations  we  have  used  2D  finite  difference  technique  which  can 
handle  very  complex  structures.  In  the  latter  ca.se  we  have  experimented  with  self-similar 
von  Karman  media,  an  uppermost  crust  low-velocity  layer  and  a  corrugated  (ID  von 
Karman  representation)  Moho  boundary.  Major  results  here  is  that  a  "scatter"  medium 
with  RMS  velocity  variations  of  c.  3-6  per  cent  are  es.sential  for  strong  coda  generation. 
Corrugations  on  major  discontinuities  are  important  in  this  context.  When  a  P-source  is 
embedded  in  a  complex  medium,  P-to-S  conversion  close  to  the  source  is  efficient.  This 
explains  why  we  see  so  much  shear  wave  energy  in  local  quarry  blast  records. 

Rg-synthetics  for  propagation  in  the  uppermost  crust  low  velocity  layer  (LVL)  are 
instructive  for  a  better  understanding  of  the  observational  recordings.  For  example,  for 
focal  depths  below  2-3  km  handly  any  Rg-waves  would  be  excited.  Thus  Rg-observations 
could  be  an  important  discriminant  for  source  classification.  However,  introducing  scatters 
in  the  LVL  tend  to  attenuate  the  Rg  besides  generating  an  abundance  of  interfering  s- 
scattering  wavelets.  In  this  way  we  can  quantitatively  explain  the  relative  elusivness  of 
Rg-pha.scs  at  least  for  stations  located  in  hilly  areas.  The  Rg-observations  used  in  analysis 
comes  from  7  arrays  on  4  continents,  and  the  estimated  LVL  shear  velocity  and  half-space 
velocity  estimates  were  remarkedly  consistent  despite  the  large  .separations  between  siting 
areas.  Average  velocities  were  2.87  kms '  amd  3.56  kms  '  respectively.  For  the  one  layer 
over  halfspace  model  used  in  inversion  of  the  dispersion  observations  the  LVL  layer 
thickness  varied  considerably  from  0.12  km  at  Yellowknife  (Canada)  to  1.60  km  at  Alice 
Springs  (Australia). 

In  section  5,  we  address  the  problem  of  event  cla.ssification  in  the  teleseismic  distance 
range  and  using  P-and  P-coda  observations  only.  Using  our  novel  discriminant  technique 
we  were  able  to  successfully  (no  mi.sclassifications)  classify  all  presumed  nuclear  tests  in 
the  E.  Kazakh  area.  However,  including  presumed  explosions  from  many  other  parts  of 
exUSSR  in  the  learning  set,  a  slight  incrca.se  of  6.5  per  cent  in  the  classifications  error 
was  observed.  This  imply  that  for  teleseismic  events  (read  explosions)  upper  mantle 
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propagation  paths  do  not  much  weaken  a  source  discriminant.  However,  for  events  with 
Mb-valucs  below  4.0  to  4..‘i  the  events  become  explosion-like  simply  due  to  p(K>r  SNR  lor 
the  coda  waves.  Since  we  used  NORKSS  observations  which  detectability  is  excellent  for 
Eurasia  we  take  the  abtwe  results  to  imply  that  Mb  4.0  -  4.5  rcprc.sent  lower  limits  for 
teleseismic  event  classification  capabilities. 

Concluding  Remarks 

During  the  contracting  period  we  had  observed  problems  related  to: 

Structural  mapping  of  the  crust/lithospherc 

Signal  detection,  event  location  and  automating  the  seismic  bulletin  production 
Wave  propagation  in  complex  media  for  realistically  synthesizing  wave  propagation 
in  the  crust/lithospherc  system. 

Extensive  experiments  with  event  classification  at  teleseismic  distances  have  also 
been  undertaken. 

Extensive  documentation  of  our  research  efforts  under  said  rcxscarch  contract,  arc  presented 
in  the  3  mandatary  annual  .science  reports  (including  this  one).  Complementary  purely 
scientific  presentations  arc  in  form  of  16  publications  listed  on  page  2  in  this  report. 

Implications  for  Further  Research 


The  need  for  a  better  understanding  of  crust/lithospherc  wave  propagation  complexities  is 
cs.scntial  for  the  design  of  llexible  event  di.scriminant.s  for  use  at  local  and  regional 
distance  ranges.  Al.so  important  is  automating  seksmograph  network  operation  so  available 
resources  can  be  focused  upon  a  few  relative  problematic  events  instead  of  wasted  on 
numerous  man-made  explosions  in  un-intcrcsting  areas. 
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Extensional  simple  shear  deformation  of  the  entire  lithosphere  has 
been  postulated  by  Wernicke  (1)  and  others  but  up  to  now 
unequivocal  seismic  evidence  in  support  of  this  hypotheses  has  been 
lacking  (2).  Here  we  report  on  well  defined  seismic  reflectors  in  the 
Skagerrak  Sea  which  are  interpreted  as  major  low  angle  shear  zones 
(The  Bambie  Fault  -  BF)  also  underlying  the  Skagerrak  graben 
segment  of  the  Oslo  Rift.  A  few  of  them  can  be  traced  almost 
continuously  from  mid-crust  through  the  lower  crust,  off-setting  Moho 
and  continuing  downwards  to  50-60  km  depth  (16  sTWT).  This  and 
similar  results  from  1730  km  of  deep  seismic  lines  in  Skagerrak 
indicate  that  the  crust  and  mantle  did  inherit  a  pronounced  structural 
fabric  from  the  earlier  Proterozoic  Sveconorwegian  orogeny. 
Reactivation  of  these  surfaces  during  the  rifting  may  explain  the 
occurence  of  simple  shear  in  the  lower  crust  below  Skagerrak. 
Conventional  rifting  scenarios  incorporating  magmatic  underplating  of 
Moho  is  not  considered  tenable  in  our  reconstruction  of  the  Oslo  Rift 
formation. 

The  Permian  Oslo  Rift  is  centrally  located  in  the  Sveconorwegian  province  of  the 
Baltic  Shield  (3,4)(  Fig  1).  The  tectonic  fabric  of  the  crust  appears  to  reflect 
alternating  E-W  compressional  and  extensional  regimes  with  the  Grenvillian- 
Sveconorwegian  orogeny  (1.2-0. 9  Ga)  as  the  most  prominent  one  (5).  For 
example,  the  Skagerrak  Sea  block  appears  to  have  been  upthrusted  and  ductily 
sheared  against  the  cratonized  Telemark  block  to  the  west  during  this  orogeny 
(6,7)  (Fig  1).  The  development  of  the  much  later  Oslo  Rift  remains  a  puzzle  despite 
much  “on  land"  geological  research  effort  (3,8).  We  have  been  able  to  adress  this 
old  problem  with  new  seismic  reflection  data  stemming  from  the  RA/  Mobile  search 
cruise  in  Skagerrak  in  1987  (9,  Fig  1).  The  entire  deep  seismic  grid  (16  s  TWT)  has 


recently  been  reprocessed  using  the  computer  facilities  at  Bullard  Labs, 
Cambridge,  UK  which  greatly  improved  the  resolution  of  the  seismic  sections. 

The  lines  having  the  most  direct  baring  on  the  rift  area  are  displayed  in  Fig. 2.  The 
sharp  reflector  in  Fig.  2a  is  outstanding,  cutting  through  the  entire  crust,  offsetting 
Moho  2-4  km  and  continuing  in  the  upper  mantle.  It  is  clearly  seen  on  both  OG-8 
and  OG-12  This  reflector  is  observable  for  over  130  km  and  is  intact  below  the 
Skagerrak  rift  segment.  We  denote  this  reflector  the  Bamble  Fault  (BF).  Its  landward 
extrapolation  coincides  with  the  PKF-fault  separating  the  Telemark  block  and 
Bamble  sector  in  Fig.  1  and  2b.  On  all  N.  Skagerrak  profiles  a  prominent  reflective 
feature  can  be  observed  originating  in  the  lower  crust  and  dipping  into  the  mantle 
(Fig  2) .  It  is  underlying  the  Skagerrak  Graben,  has  a  highly  reflective  character  and 
its  upper  surface'  in  the  mantle  is  mapped  out  in  Fig  3a.  We  interpreted  it  as  the 
remains  of  the  underthrusted  crustal  segment  of  the  Telemark  craton.  It  is  clearly 
shown  on  both  intersecting  profiles  OG-8  and  OG-2  (Fig.  2b,c).  The  structural  style 
of  this  segment  is  similar  to  the  crust/mantle  imbrications  that  are  observed  below 
young  orogenic  zones  like  the  Pyrenees  and  the  Alps  (10).  The  Moho  disruption  in 
the  southwest  on  profile  OG-2  (Fig.  2c),  coincides  with  the  Fjerritslev  fault  rone 
(FFZ;  Fig.  1). 

Pronounced  crustal  thinning  coincides  with  the  Oslo  Rift  axis  except  in  the  south 
where  the  Skagerrak  graben  axis  is  shifted  s»ightly  westward  (Fig.  3b). 

Dextral  movements  (ca.  3  km)  are  reported  along  the  PKF-fault,  bounding  the 
Bamble  sector,  part  of  which  took  place  in  Permian  times  (4).  Palaeozoic  sediments 
have  been  reported  in  the  entire  Oslo  Rift  (11).  presumed  deposited  in  an  older 
pre-rift  basin.  Outside  the  Skagerrak  Graben  only  post-rift  sedimentary  strata  are 
found  (Fig.  3c).  In  our  tectonic  scenario  the  Skagerrak  block  was,  during  the 
Sveocnorwegian  orogeny,  fragmented  with  eastward  parts  being  slightly 
overthrusted  the  westward  ones  (Fig.  2b).  this  is  similar  to  the  observations  of  the 
crustal  shortening  allong  the  contemporary  Grenvillian  front  in  eastern  Canada 
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(12.13). 


In  summary,  the  general  Oslo  Rift  area  was  in  Precambrian  times  subject  to  both 
compressional  arKl  extensional  deformations  thus  weakening  the  lithosphere.  The 
observed  thinned  crust  of  this  area  (Fig  3a)  cannot  alone  be  ascribed  to  the 
moderate  Permian  extension.  Dyke  and  granite  intrusions  of  1000  Ma  and  younger 
in  the  Skagerrak  coastal  areas  (8.14)  also  imply  a  lithosphere  weakening. 

During  the  early  Permian  extension  the  landward  and  seaward  parts  of  the  Oslo 
Rift  apparently  responded  very  differently.  According  to  the  numerical  extension 
models  of  Dunbar  and  Sawyer  (15)  the  continental  lithosphere  is  analogue  to  a 
composite  material  of  alternating  high-  and  low-strength  layers.  A  weak  lower  crust 
would  act  as  a  detachment  zone  relatively  to  a  strong  upper  crust  and  a  strong 
upper  lithosphere.  For  a  region  being  subjected  to  extensional  forces  ,  initial  rifting 
would  be  located  in  areas  having  pre-existing  weaknesses.  In  the  northern 
onshore  segment  of  the  Oslo  Rift  vertical  faulting  dominated,  cracking  the 
weakened  lithosphere  with  subsequent  influx  of  basaltic  magmas  in  the  initial 
stages  of  IQ-20  Ma  (3.8).  Stretching  was  very  modest  (8  approx.  1.1)  and  no  pre-rift 
dooming  (16).  In  the  N. Skagerrak  simple  shear  deformation  apears  to  have  taken 
place  along  the  Bamble  Fault  (Fig.2b)  and  with  the  crust  /  lithosphere  acting  as  a 
uniform  block.  Extensional  reactivation  of  old  Sveconorwegian  compressional 
shear  zones  in  the  upper  and  lower  crust  would  explain  the  favoured  motion  along 
a  localized  faultplane  through  the  lower  crust  (15).  The  lower  crust  is  othenvise  by 
rheological  and  temperature  considerations  presumed  to  deform  by  distributed 
plastic  stretching  and  thus  not  to  accomodate  simple  shear  deformation.  The 
survival  of  the  Precambrian  “Telemark”  underthrust  beneath  the  Skagerrak  Graben 
is  taken  in  support  of  this  view.  Further  south  deformations  may  have  taken  place 
by  pure  shear  along  detachment  zones  in  the  lower  crust  and  also  in  the  lower 
lithosphere  (2,17)  For  deep  detachment  zones  in  the  lowermost  crust,  the 
corresponding  crustal  thinning  would  be  small  arxl  shifted  relative  to  the  graben 
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axis  as  actually  observed  (Fig  3b.c).  Likewise,  geometry  of  the  post  rift  sedimentary 
basin  is  also  observed  to  be  shifted  (Fig.  3c).  Our  Oslo  Rift  models  are  cartoonized 
in  Fig.  4. 

Dominance  of  simple  shear  deformation  would  give  rise  to  rotation,  for  which  there, 
as  mentioned,  is  evidence  along  the  Bamble  sector  (Fig.  3b).  The  Skagerrak  block 
is  detached  along  the  FFZ  (18)  (Fig.1  and  2c)  where  Permian  magmatic  activity  has 
been  reported  (29,20).  Finally,  substantial  parts  of  the  lithosphere  cannot  have 
been  distorted  too  much  as  deep  reflectors  down  to  around  100  km  were  reported 
for  the  OG-13  line  in  Fig.1  (21) 

We  have  outlined  the  major  trends  in  the  Oslo  Rift  formation  while  on  a  smaller 
scale  length  secondary  deformation  features  in  Skagerrak  have  also  been  reported 
(22).  We  find  no  evidence  of  Moho  magmatic  underplating  in  Skagerrak  a 
phenomenon  often  postulated  for  rift  zones.  The  outstanding  observational  features 
are  the  tracing  of  the  BF  down  to  at  least  60  km  and  the  remains  of  the 
underthrusted  Telemark  craton  beneath  the  graben.  We  conclude  that  our 
observations  strongly  indicate  that  simple  shear  deformation  of  the  lithosphere 
occured  below  the  N.  Skagerrak  and  that  the  mode  of  lithosphere  deformation  is 
influenced  by  pre-existing  structural  fabric. 
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6°  8°  10“  12“ 


6“  8“  10“  12“ 


Fig  la.  Profiling  grid  -  RA/  Mobil  Search  cruise  in  Skagerrak,  1987. 
Aquisition  parameters:7000  cu  Inch  airgun  array,  50  m  shotpoint  interval, 
4500  m  long  streamer  with  50  m  group  spacing  and  4  ms  samplerate. 
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Fig  1b  Tectonic  map  of  the  southwestern  part  of  the  Baltic  Shield  based 
on  works  of  Gaal  and  Gorbatschev  (5)and  Kinck  et  al  (14)  AG:  Akershus 
Graben,  BF;  Borgum  Fault,  BZ  Bamble  Zone,  CDF:  Caledonian  Deformation 
Front.  FBZ  Fennoscandian  Border  Zone.  FFZ:  Fierritslev  Fault  Zone,  LF 
Loten  Fault  Zone,  MZ  Mylonite  Zone,  MAF  Meheia-Adal  Fault  Zone,  NDB 
Norwegian-Danish  Basin,  OF:  Oppland  Fault  Zone,  PKF:  Porsgrunn- 
Kristiansand  Fault  Zone.PZ:  Protogine  Zone,  SG:  Skagerrak  Graben,  SNF: 
Sveconorwegian  Front,  TSIB:  Trans-Scandinavian  Ignious  Belt.  VG:  Vestfold 
Graben,  WGR  Western  Gneis  Region,  OF  Oymark  Fault 
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neament ,  marked  with  an  arrow,  dipping  through  the  crust  and  mantle  from  partly  masked  by  strong  water  bottom  multiples,  but  is  clearly  identified  on 

the  gravity  profile 


r.ClMOn^S  71WIL-A.VM-OAM 


2c)  Interpreted  migrated  section  ot  profile  OG-2  (Fig  1).  Note  the  4-5  km  Skagerrak  graben  and  Moho  is  marked.,  in  the  notheast  the  underthrusted 

offset  in  Moho  directly  below  the  Fjerritslev  Fault  Zone  (FFZ)  indicating  segment  of  the  Telemark  craton  is  marked  with  (A),  as  on  Fig  2  a, b. 

vertical  faulting  cutting  the  entire  crust..  Base  Triassic  unconformity  (BTU),  ,  , 


Fig  3 


3a)  Contours  of  the  upper  surface  of  the  underthrusted  segment  of  the 
“Telemark  Craton”  in  the  mantle  shown  in  Fig.  2  marked  (A). 

3b)  Thickness  of  crystalline  crust  in  Scandinavia  (6).  Contour  interval  2 
km,  the  Skagerrak  Graben  contours  are  indicated. 

3d)  Interpreted  migrated  section  of  line  OG-7;  note  that  the  post-rift 
sediments  are  here  shifted  eastward  relative  to  the  graben  axis. 

I  3 


Post-rift  Pre-rift 


Northern  rift  segment.  Central  rift  segment,  Southern  rift  segment, 

landward  N.  Skagerrak  S.  Skagerrak 


West  East  West  East  West  East 


Schematic  view  of  'he  Oslo  Rift  evolution  for  its  northern,  central  and 
southern  graben  segments.  Dimensions  not  to  scale.  The  upper  figures 
reflect  the  pre-rift  stages  and  the  lower  reflect  the  post  rift  stages.  Shading 
indicate  upper/lower  crust  and  likewise  for  the  mantle  lithosphere 
Northern  segment.  (Fig  4  a):  Paleozoic  pre-rift  sedimentary  basin  implying 
thinned  crust  and  weakened  lithosphere  (Fig  4b);  Rifting  with  graben 
formation  and  extensive  magmatism.  Moderate  stretching  Pre-rift  sediments 
only  preserved  within  the  graben  due  to  base  Tertiary  erosion  Central 
segment,  (Fig  4c)  The  pre-rift  crustal  fabric  dominated  by  the 
Svecorxjrwegian  crustal  shortening  (Fig  4d)  Simple  shear  deformation  of 
the  entire  lithosphere  No  magmatism  observed,  less  dominant  graben 
formation  Post-rift  sedimentation  and  maximal  crustal  thinning  slightly 
shifted  to  the  east  Southern  segment,  (Fig  4e)  Weakened  Precambrian 
crustal  imprint  relative  to  the  central  segment  (Fig  2  a,b)  observed.  (Fig  4f) 
Deformation  by  pure  shear  in  the  lower  crust  with  detachment  slightly  offset 
to  the  east,  explaining  asymmetry  between  graben  axis  and  general  crustal 
thinning  and  post-rift  sedimentation  Well  developed  graben  formation  with 
no  dear  indications  of  magmatism  Magmatism  is  however  reported  further 
to  the  south  along  the  FFZ  (Fig  1  and  Fig  2c) 


Fig  4 
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Abstract 

In  general,  the  continental  upper  crystalline  crust  appears  almost 
nonreflective  on  many  deep  seismic  reflection  profiling  sections.  Only  major 
faults  and  lithological  boundaries  can  be  localised,  while  smaller  scale  faults 
are  usually  not  identified.  A  deep  seismic  survey  covering  the  Skagerrak  sea 
was  undertaken  in  1987  by  the  RA/  Mobil  Search  -  depth  range  0  -16  s  TWT 
Although  the  north-eastern  part  of  the  Skagerrak  has  only  thin  or  no 
sediment  coverage,  the  data  in  the  depth  range  0  -  6  s  TWT  were  still 
reprocessed  to  commercial  exploration  standard.  Close  inspection  of  these 
high  quality  sections  provides  evidence  of  seismic  images  of  an  abundance 
of  moderately  steep  (20-50  deg)  faults  that  can  be  traced  to  depths  of  10-13. 
Fault  geometry  can  be  inferred  at  profile  intersections,  and  coincides  with 
the  general  tectonic  fabric  observed  on  land  and  beneath  the  Skagerrak 
Sea.  These  areas  were  deformed  during  the  Proterozoic  Sveconorwegian 
orgenyin  and  the  formation  of  the  Permian  Oslo  Rift.  Due  to  the  post-Permian 
flattening  of  the  Skagerrak  Sea  proper,  the  relative  age  and  movements  of 
the  basement  faults  could  not  be  given.  This  study  shows  that  faults  in  the 
upper  crystalline  crust  are  mappable  by  seismic  means  under  favourable 
conditions.  Such  measurements  have  the  potential  of  resolving  such 
problems  as  block  geometries  with  depth  and  recent  seismotectonic 
movements 


16 


Introduction: 


Deep  seismic  reflection  profiles  from  the  continents  often  show  a  reflective 
lower  crust  and  Moho.  In  comparision  .the  upper  crystalline  basement  (down 
to  10-20  km  depth)  usually  lacks  coherent  reflections  and  in  turn  is 
dominated  by  scattered  energy.  This  statement  does  not  exclude  presence 
of  upper  crustal  reflectors  which  typically  are  observed  in  areas  of  major 
crustal  deformation  such  as  the  Grenvillian  province  in  eastern  Canada. 
(Green  et.  al.,  1989)  Similarly,  large  scale  faults  in  the  upper  crystalline  crust 
are  easily  identified  in  seismic  reflection  profiles  (Klemperer  and  Hurich, 

1990,  Bois,  1991,  Brewer  et.  al.  1980)  This  paper  addresses  seismic 
identification  of  relative  small-scale  faults  in  basement  rocks,  which  may  be 
related  to  those  observed  in  exposed  rocky  terrain. 

In  1987  a  deep  seismic  survey  was  conducted  in  the  Skagerrak  Sea 
by  the  R/V  Mobil  Search  (Fig.1,  acquisition  details  in  Table  1).  The  original 
onboard  processed  seismic  sections  revealed  numerous  lower  crustal  and 
upper  mantle  features  (e  g.  see  Husebye  et  al.l988,  Larsson  and  Husebye 

1991,  Lie  et  al.,  1990  and  Pedersen  et  al.  1990).  The  original  onboard 
stacked  data  (0-16  s  TWT)  were  subjected  to  a  new  extensive  post-stack 
processing  sequence  including  FK-migration  (details  in  Table  3)  as  a  part  of 
a  cooperative  study  with  the  British  Institutions  Reflection  Profiling  Syndicate 
(BIRPS).  To  utilize  the  full  information  potential  of  the  collected  seismic  data 
the  upper  6  sec.  of  the  sections,  limited  to  0-6  sTWT,  were  reprocessed  to 
commercial  exploration  standard  including  DMO  and  post-stack  migration 
(processing  details  in  Table  2). 

The  enhanced  quality  of  the  reprocessed  seismic  sections  has 
enabled  us  to  investigate  some  fundamental  aspects  of  marine  seismic 
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profiling  resolution  such  as  the  extent  to  which  faults  in  the  crystalline 
basement  are  seismically  mappable.  Commercial  exploration  seismic 
surveys  have  rarely  been  attempted  in  areas  with  marginal  sediment  cover, 
as  in  the  northern  Skagerrak  Sea.  On  the  following  high-quality  profiles  , 
seismic  images  of  subtle  but  consistent  faulting  in  the  upper  crystalline  crust 
are  described. 

The  Skaoerrak  tectonic  setting 

The  Skagerrak  Sea  is  considered  to  be  a  continental  platform  area  between 
the  Sveconorwegian  terrains  of  the  stable  Baltic  Shield  to  the  northwest  and 
the  mobile  belts  of  north-western  Europe  to  the  southwest.  The  separation  of 
the  two  tectonic  units  is  along  the  Fennoscandian  Border  Zone  (FBZ)  as 
shown  in  Fig.1.  The  landward  Sveconorwegian  rocks  are  of  Proterozoic  age 
(1.7  -1.4  Ga)  and  the  youngest  are  the  Permian  volcanics  of  the  Oslo  Rift 
(Gaal  and  Gorbatchev,  1987,  Bjorlykke  et.  al..  1990).  The  Sveconorwegian- 
Grenvlllian  orogeny  (1.25-0.85  Ga)  affected  the  southern  part  of  Norway  and 
Sweden  with  extensional  and  compressional  forces  mainly  acting  in  the 
E/W-direction  (Falkum,  1985).  The  onshore  Bamble  zone,  a  high  grade 
metamorphic  terrain,  is  hypothesized  to  have  been  uplifted  an  estimated  15- 
20  km  relative  to  the  westward  Telemark  block  along  the  Porsgrunn- 
Kristiansand  fault  zone  during  the  above  orogeny  (Falkum  and  Starmer, 
pers.  comm  ).  On  the  deep  seismic  profiles  major  structural  elements  can  be 
identified  in  the  crust  and  mantle  under  the  N.  Skagerrak  Sea  which 
correlate  well  with  the  geometry  of  the  tectonic  boundaries  mapped  on  land 
to  the  northwest  (Lie  and  Husebye,  1992,  Lie  et  al.  1992b). 

There  is  little  evidence  of  tectonic  activity  in  the  late  Precambrian  ,  but 
probably  the  general  area  was  of  low  relief  (Ramberg  and  Spjeldnaes, 
1978)  The  Permian  graben  formation  of  the  Oslo  Rift  was  proceeded  by  the 
formation  of  narrow  sedimentary  troughs  in  Cambro-Silurian  times 
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(Bjorlykke,  1983).  Extensive  volcanism  within  the  landward  part  of  the  Oslo 
Rift  ( the  Vestfold  and  the  Akershus  grabens  )  is  obvious  but  the  rift  appears 
to  tack  a  significant  volume  of  volcanics  in  its  seaward  extension,  the 
Skagerrak  Graben  (Fig.  1)(  Lie  et  al.,  1992b).  Erosion  in  Late  Permian  /  Early 
Triassic  times  has  affected  the  entire  Skagerrak  Sea  proper  with  Paleozoic 
sediments  preserved  only  within  the  Oslo  Rift  system. 

The  major  faults  in  the  Skagerrak  Sea  are  shown  in  Fig.  1 .  The  general 
tectonic  trend  direction  is  N/S  (NNE  and  NNW  in  the  coastal  areas  to  the 
east  an  west,  respectively)  (Ramberg  and  Gabrielsen,  1977).  The  major 
onshore  faults  appear  to  have  formed  prior  or  during  the  middle  Proterozoic 
Sveconorwegian  orogeny  (Starmer,  1985)  but  some  have  clearly  been 
reactivated  at  much  later  times,  especially  during  the  Permian  rift  phase 
(Spjeldnaes,  pers.  comm  ).  The  landward  faults  are  well  established  as 
major  structures  on  geological  grounds,  but  fault  geometries  at  depths 
greater  than  a  few  kilometres  are  uncertain.  In  the  northeast  Skagerrak 
Sea,  faults  have  been  mapped  essentially  on  the  basis  of  bathymetric  and 
shallow  seismic  (sparker)  data  (Floden,  1973).  In  the  southern  Skagerrak 
Sea  the  OG-profiles  and  other  exploration  seismic  profiles  have  been  used 
to  map  the  Skagerrak  Graben  (  e  g.  see  Husebye  et  al.,  1988,  Lie  and 
Husebye,  1992,  and  Kinck  et.  al.,  1991  among  others). 

Seismic  identification  of  faultolanes  in  the  crystalline  crust. 

The  deep  seismic  profiles  from  Skagerrak  (Table  2)  provide  clear  crustal 
cross-sections  as  shown  in  Fig  2.  The  bounding  faults  of  the  Skagerrak 
Graben  (Fig. 2)  form  strong,  dipping  reflections  that  are  easily  recognized 
seismically,  similar  features  documented  from  many  other  places  in  the 
world.(e.g.  see  Leven  et  al.,  1990).  However,  the  smaller  scale  faults  that 
outline  the  basement  blocks  on  the  surface  of  rocky  terrain  are  not  easily 
identified  (Fig.  2,  line  OG-13).The  seismic  signature  of  the  upper  crystalline 
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crust  (say  above  5  sTWT)  is  often  dominated  by  scattered  energy  and 
somewhat  chaothic  reflections.  Close  inspection  of  the  migrated  sections, 
reveals  dipping,  weak  seismic  features  in  the  upper  crystalline  crust  (fig.  2- 
5).  These  examples  are  acquired  from  areas  where  the  sediment  cover  is 
less  than  1-2  km.  We  interpreted  these  features  as  moderately  dipping  faults. 
Supportive  evidence  for  this  interpretation  is  that  these  presumed  faults  can 
be  identified  jointly  at  both  sections  at  intersecting  profiles  as  demonstrated 
in  Fig.  4.  Their  strike  generally  coincides  with  the  main  tectonic  lineaments  of 
the  surrounding  onshore  area,  and  also  with  the  inferred  fault  trends  from 
early  bathymetry  and  seismic  sparker  data  (Floden,  1973).  The  faults  in  the 
crystalline  basement  are  in  places  observed  to  coincide  with  offsets  in  the 
overlaying  sedimentary  cover  (Figs.  4  and  5)  The  fault  planes  can  be 
mapped  down  to  4-5  s  on  the  time  sections  equivalent  to  approximately  10- 
13  km  depths  and  their  appearent  dip  angles  are  in  the  range  20-50  deg. 

(Fig.  6). 

Steep  fault  planes  are  generally  attenuated  in  a  stack  which 
emphasizes  horizontal  reflectors  In  sedimentary  strata  such  faults  are 
indirectly  identifyed  through  abrupt  offsets  in  the  stratification,  and  not  as  a 
direct  image  of  the  fault  plane  itself.  Similarly,  in  the  chaotic  pattern  of  the 
upper  crust,  steep  faults  are  not  only  imaged  as  reflections  from  the  fault 
plane,  but  indirectly  as  aligned  terminations  of  short  diffractive  “events” 
against  the  fault  The  fault  planes  are  barely  visible,  if  only  viewed  as  offsets 
in  the  uncorrelated  "noise”  at  one  “time-sample  level",  but  as  alignments 
across  the  section  they  are  more  easily  identified.  The  fault  signature 
discussed  above  is  subtle  and  can  easily  be  overlooked  especially  when  the 
quality  of  the  seismic  sections  is  not  optimum.  In  either  version  of  the 
reprocessed  sections  (  Table  2  and  3).  the  post-stack  migration  of  the  data 
has  had  a  major  piositive  effect  on  bringing  out  fault  plane  images,  although 
some  also  appear  on  the  unmigrated  sections  Reflections  from  minor  faults 
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in  crystalline  basement  have  also  been  reported  from  Sweden  where  the 
observations  later  were  confirmed  by  drilling  (Juhlin,  1990).  The  bore-hole 
log  information  showed  that,  the  thin  layers  of  fractured  granite  gave 
adequate  impedance  contrasts  to  provide  observable  reflections. 

Many  of  the  above  described  dipping  crustal  features  have  weak 
signatures,  so  alternative  explanations  for  these  observations  have  been 
considered.  The  most  obvious  candidates  are  diffractions  and  reflected 
refractions  (Young  et.  al..l990).  We  tested  this  hypothesis  by  using  models 
of  zero-offset  diffraction  hyperbolas  based  upon  representative  velocity 
functions  (Fig.  7)  These  in  turn  were  compared  to  the  interfered  faults  but 
did  not  coincide.  Out-of-plane  diffractions  are  more  difficult  to  verify  or 
exclude  as  they  may  appear  within  a  much  wider  dip  range.  Potential  reflect- 
refraction  events  were  also  modelled  (Fig  7)  but  did  not  provide  a  likely 
alternative  to  the  fault  interpretations  in  Figs  2-5.  To  summarize,  we  claim 
that  the  above  seismic  sections  /  structural  examples  are  probably  not 
attributable  to  diffractions  or  processing  artefacts.  In  essence,  we  consider  it 
feasible  to  identify  by  seismic  means  faults  in  the  upper  crystalline  crust 
given  good  profiling  data  in  areas  with  sparse  sedimentary  coverage 

Basement  Fault  Mapping 

Crystalline  basement  faults  are  most  commonly  observed  in  the  northern 
Skagerrak  Sea  where  the  sediment  cover  is  less  than  1-2  km.  Hence  our 
mapping  efforts  were  limited  to  the  profiles  OG-8,  OG-9,  OG-12,  OG13  and 
the  northern  segments  of  OG-2  and  OG-3.  The  results  are  shown  in  Fig  6  in 
the  form  of  time-to-depth  converted  line  drawings  using  a  representative 
velocity-depth  function.  The  fault  density  is  relatively  high  with  a  separation 
distances  amounting  to  a  few  kilometres.  Fault  intersections  are  frequently 
observed,  albeit  the  dominant  trend  is  towards  the  Skagerrak  Graben. 
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Implicit  in  the  later  statement  is  that  the  identified  faults  are  roughly  parallel 
to  this  dominant  structural  feature.  Evidence  for  this  in  terms  of  determined 
strike  and  dip  of  fault  planes  is  only  available  around  line-crossings  like  the 
OG-13/OG-8  as  demonstrated  in  Fig  4  Similar  strike  directions  were  found 
at  all  other  line  crossings  in  the  northeastern  Skagerrak 
Sea 

The  relative  movements  along  the  mapped  basement  faults  are 
difficult  to  assess  because  the  late  Permian  /  early  Triassic  erosional  event 
post-dates  the  latest  tectonic  deformation  episode.  This  erosional  event, 
following  the  Permian  rifting,  removed  the  Paleozoic  sediment  cover  in  the 
Skagerrak  Sea  proper  and  flattened’  the  basement  topography  except  in 
the  Oslo  Rift  (Fig.  2).  The  overlying  thin  Mesozoic  and  Cenozoic  sediment 
cover  in  the  northern  Skagerrak  Sea  is  relatively  undeformed.  However,  in 
the  southern  Skagerrak  Sea  which  has  a  thick  sedimentary  cover,  we  were 
unable  to  identify  basement  faults.  An  exception  here  is  the  low  angle 
bounding  faults  of  the  Skagerrak  Graben  (Fig  2)  where  Cambro-Silurian  and 
syn-rift  sediment  thickness  may  exceed  3-4  km.  Note,  in  the  latter  case  these 
faults  where  observed  even  in  the  original  stacked  sections. 

Discussion 

In  the  preceeding  sections  we  have,  in  our  opinion,  provided  evidence 
moderately  dipping  faults  can  be  imaged  Beneath  the  Skagerrak  Sea  one 
precondition  for  this  is  that  the  sediment  cover  does  not  exceed  a  few 
kilometres  which  would  lower  the  signal-to-noise  ratio  The  data  resolution 
can  be  adversely  effected  by  scattering  effects  in  areas  with  a  hard  uneven 
sea  floor  with  no  sediments  (Lamer  et.  al.,  1983,  Lie  et.  al.,  1990).  Our 
interpretation  of  these  fault  signatures  is  supported  by  the  fact  that  in  places 
they  do  coincide  with  offsets  in  the  sediments  and  sea  floor  Fault  planes  can 
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also  be  mapped  on  crossing  profiles  and  their  infered  strike  and  dip 
coincides  well  with  the  overall  tectonic  trend  of  the  surrounding  area. 

Floden  (1973)  and  others  have  on  basis  of  sparker  profiles  and 
bathymetric  mappings  “identified”  a  few  prominent  faults  in  the  northeastern 
Skagerrak  Sea  (Fig  i)  These  hyothesized  faults  are  not  uniquely 
recognized  in  our  results  (Fig  6)  This  may  be  due  to  the  relatively  long 
source-near  trace  offset  on  the  Mobil  Search  lines  which  in  shallow  waters 
yields  poor  resolution  in  the  uppermost  part  of  the  section.  This  in  turn 
makes  correlation  to  detailed  bathymetric  data  difficult 

The  northern  Skagerrak  sea  results  may  be  used  to  resolve  tectonic 
evolution  problems.  For  example,  most  of  the  graben  cross-sections,  with 
exception  of  profile  OG-7  (Fig  2),  show  asymmetric  half-grabens  with 
prominent  bounding  faults  The  main  strike  of  these  faults  is  mostly  to  NNE 
which  roughly  coincides  with  the  infered  strike  directions  of  the  mapped 
basement  faults.  On  the  basis  of  our  findings  we  conclude  that  the 
extensional  deformation  during  the  Permian  rift  phase  has  not  been  limited 
to  the  narrow  Skagerrak  Graben  per  se.  but  probably  was  absorbed  along 
the  many  basement  faults  mapped  in  the  entire  NE  Skagerrak.  This  has 
importance  for  properly  estimating  the  crustal  extension  during  the  rifting 
(Pedersen  et  al.,  1991)  Note,  the  Permian  Oslo  Rift  is  located  in  an 
inherited'  Sveconorwegian  structure  implying  that  basement  faults  could 
either  be  the  same  ageas  the  rift  or  most  likely  considerably  older  (Ramberg 
and  Spjeldnaes,  1978,  Falkum,  1985). 

Finally,  we  want  to  address  the  so-called  seismotectonic  concepts 
popular  in  earthquake  hazard  analysis  (eg  see  Gregersen  and  Basham, 
1989).  These  kinds  of  problems  are  of  interest  in  the  Skagerrak  Sea  as  the 
largest  Fennoscandian  earthquake  in  modern  times  (Husebye  et  al.,  1978) 
took  place  in  the  outer  Oslo  Fjord  in  1904.  Occasionaly  we  observe  that 
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basement  faults  offset  the  present  day  Mesozoic/Cenozoic  sediment  cover  in 
the  northern  Skagerrak  Sea  implying  that  the  area  has  not  remained 
tectonically  quiscent.  In  a  seismotectonic  context  fault  candidates  are 
numerous;  the  graben  bounding  faults  .  one  of  the  many  faults  outside  the 
graben  or  the  prominent  low  angle  Bamble  shear  zone  extending  from  the 
Bamble  sector  towards  the  the  Swedish  west  coast  (Lie  et  al.,  1992a). 
Although  the  epicentre  of  the  1904  earthquake  is  not  well  constrained,  the 
point  is  that  under  favourable  conditions  if  should  be  possible  fo  map 
seismically  active  faults  with  reflection  seismic  profiling.  High  resolution 
profiling  systems  are  needed  to  resolve  fault  offset  in  the  youngest 
sedimentary  strata. 

Concluding  remarks 

We  have  reported  on  observations  of  moderately  dipping  faults  within  the 
upper  crystalline  crust  on  reflection  profiles.  Such  faulfs  are  nof  generally 
recognized  seismically  because  full  exploration  standard  marine  seismic 
surveys  are  not  often  undertaken  in  areas  with  little  sedimentary  cover.  The 
subtle  fault  signatures  documented  here  are  therefore  easily  overlooked  or 
dismissed  as  artifacts.  As  conventional  seismic  profiling  strategy  and  data 
processing  are  not  focused  on  non-horizonfal  structural  features  our 
interpretatoins  may  be  somwhaf  uncertain  because  of  this  basic  limitation.  In 
a  later  paper  we  will  address  this  problem  in  terms  of  generating  synthetic 
reflection  seismograms  and  using  commercial  software  for  stacking  and 
migration  The  usefulness  of  such  approaches  have  been  demonstrated  by 
Blundell  and  Raynaud  (1986)  ,  Raynaud  (1988)  and  Cao  ef  al,  (1990) 
among  others 
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RECORDING  PARAMETERS, 

R/V  Mobil  Search  survey  Skagerrak  1987 


RECORDED  BY: 

MOBIL  ESP  SERVICES  INC 

INSTRUMENT  TYPE: 

Tl  DFS  5.5 

RECORD  LENGTH: 

16  000  MSEC 

SAMPLE  RATE: 

4  MSEC 

TRACES/RECORO. 

180 

FIELD  FILTER: 

LOWCUT  5  3  HZ 

HIGHCLfT  90  HZ 

SLOPE: 

LO  - 18  DlVOCT 

HIGH  -  72  (ju/OCT 

HYDROPHONES  PR  GROUP 

32 

GROUP  INTERVAL 

25  M 

NUMBER  OF  GROUPS: 

180 

ACTIVE  STREAMER  LENGTH 

4499  M 

STREAMER  DEPTH 

15  M 

OFFSET. 

262  M 

1  TUNED  WIDE  AIR  GUN  ARRAY  I 

ARRAY  DIMENTIONS: 

59  M  LONG  X  75  M  WIDE 

NO.  OF  SUB-ARRAYS 

4  X  (12  GUNS  OVER  59  M) 

SUB-ARRAY  SPACING: 

25  M  BETWEEM  STRINGS 

PRESSURE/CAPASITY 

1800  PSI/7200  CU  IN 

TOWING  DEPT: 

10M 

SHOTPOINT  INTERVAL: 

50  M 

CDP  SPACING: 

12.5  M 

FOLD 

45 

Table  1.  Recording  parameters  during  the  R/V  Mobil  Search  seismic 


profiling  cruise  in  Skagerrak  Sea,  1987. 
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PROCESSING  SEQUENCE,  R/V  Mobil  Search  survey 
Skagcrrak  1987,  Reprocessed  0-16  s  TWT,  1990/1991 


(laS;.  ONDOAHD  »/V  MOim  StAHCH) 

1  INSTHUMt  NT  DUniASl  (  II  ILH 

2  TRACE  MUTt 

3  GAIN  CONTROL 

4  TRACE  COWPOSiTt' .  5  TO  1.  SEVEN  WINDOWS 

5  GAPED  DECONVOLUTION 

GATE  200  6000  s  1  ENGTH  248  ms 
GAP  =  32  ms 

GATE  6000  13000  s  LENGTH  352  ms 
GAP  ^  64  ms 

6  CDP  GATHER 

7  NMO 

Velocity  .liulysts  every  3  Em 

8  TRACE  MU  1 1 

9  M(  AN  stack  43EOII) 

with  spike  reicclioii 


(1990'1091  I’o'.t  Slack  leprocessiiif).  ETlTtPS,  Cambridge,  UK; 


1  resample  IE)  ms 

2  TRACE  SUMATION 

2  adjacofil  traces 

3  GAIN  BALANCING 

varying  design  and  application 

4  GAPE D  DECONVOLUTION 

GATE  1400  5000  s  LENGTH  248  ms 
GAP  n  36  ms 

GATE  6000  13000  s  1  ENGTH  fcOO  ms 
GAP  =  85  ms 

5  FK  OIP-FIETER 

polygone  lejecling  steeply  dipping  events 

6  AGC  SCALING 

window.  250  ms  gale  50  ■  4C00  ms  pci  20 

window.  1000  mr,  gale  50  6000  ms  pet  30 

window.  5000  ms  gale  50  -  16000  ms  pel:  40 

7  EK  migration 


voioc*iy  !■ 

ffio  hini.hon 

H  fiANf  if'Ar.-.  f 

II  1)  It 

1  < 

■ 

h.  sl..|„ 

1 

I.lt.'r.it; 

Ih.-l 

1  H 

.'4  HI 

9  l)l:aT  AY 

bias  (.7  g.iin  2 


Table  2.  Processing  sequence, for  the  reprocessing  of  the  0-6  s  TWT  seismic 
data  from  RA/  Mobil  Search  cruise  in  Skagerrak  Sea,  1987. 
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PROCESSING  SEQUENCE,  R'V  Mobil  Search  survey 
Skagerrak  1987,  Reprocessed  0-6  s  TWT,  1989 


1 .  FK-FILTERING  ON  SHOT  RECORDS 

Dips  rejected  *14  to  -  5  ms/Irace 
passed  :  ♦  8  to  -1  ms/trace 

2  SPHERICAL  SPREADING  COMPENSATION 

3  DRS 

Gaped  deconvolution 

Gap  32  ms  Filter  length:  240  +  32  ms 

4  COP  SORT 

4S  told  12  5  m  CDP  distance 

5  VELOCITY  ANALYSIS 

6  SCALING 

AGC  scaling  pre  OMO.  Window  ;  400ms 

7  OIPMOVEOUT 

NMO  with  stacking  velocities 
correction  on  45  ollsel  planes 

8  01  SCSCALE 

Removal  ol  previous  400ms  AGC  alter  unNMO 

9  NMO  WITH  STACKING  VELOCITIES 
to  INNER  AND  OUTER  TRACE  MUTE 
It  SCALING 

AGC  scaling  pre  slack,  window;  250  ms 

12  STACK 

45  told  stack.  t2  5  m  CDP  distance 

13  MIGRATION 

F inile  dilterence  lime  migration 
migration  layer  thickness  40  ms 

14  DIP  FILTER 

Rejeaed  ♦  3  5  to  -3  5  ms/trace 
Passed  +  2  to  -2  ms/trace 

15  DAS 

Gaped  deconvolution 

gap:  32  ms  Filter  length:  240  ♦  32  ms 

16  LINEAR  GAIN 

Varying  design  and  application 

17  TIME  VARIANT  FILTER 


TIME 

LC 

SLOPE 

HC 

SLOPE 

(msl 

Ih7] 

Idb/ocl) 

|0/| 

(db/oct) 

0  500 

18 

72 

70 

4P 

750  1000 

1? 

48 

70 

4ft 

4000  4200 

H 

48 

40 

7.? 

5000  (.000 

n 

48 

30 

7i» 

1H  :  f'Al  tNG 

50  fobusi  AOC 
19  (DISPLAY 


Table  3.  Processing  sequence, for  the  reprocessing  of  the  0-16  s  TWT 
seismic  data  from  FW  Mobil  Search  cruise  in  Skagerrak  Sea,  1987. 
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Fig  1a.  Profiling  map,  RA/  Mobil  Search  cruise  in  Skagerrak,  1987. 
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Fig  1  b.  Tectonic  map  of  the  southwestern  part  ot  the  Baltic  Shield  based  on 
works  of  Gaal  and  Gorbatschev  (1987),  EUGENO-S  Working  Group  (1988), 
Lind  (1982)  and  Floden  (1973).  AG:  Akershus  Graben,  BF;  Borgum  Fault, 
BZ:  Bamble  Zone,  CDF:  Caledonian  Deformation  Front,  FBZ: 
Fennoscandian  Border  Zone.  FFZ:  Fierritslev  Fault  Zone,  LF:  Loten  Fault 
Zone,  MZ  :  Mylonite  Zone.  MAF:  Meheia-Adal  Fault  Zone.  NDB;  Nonwegian- 
Danish  Basin.  OF:  Oppland  Fault  Zone.  PKF:  Porsgrunn-Kristiansand  Fault 
Zone.PZ:  Protogine  Zone.  SG;  Skagerrak  Graben.  SNF:  Sveconorwegian 
Front.  TSIB:  Trans-Scandinavian  Ignious  Belt.  VG:  Vestfold  Graben.  WGR: 
Western  Gneis  Region.  0F;  Oymark  Fault.  Note,  the  Oslo  Rift  comprises  the 
Skagerrak.  Vestfold  and  Akershus  Grabens. 
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TWO-WAY-TIME  SECONDS 
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Fig  3  Detail  A  of  profile  OG-13  shown  in  fig  2.  The  most  prominent  fault',  n  also  marked  with  an  arrow  at  i  5  s  and  horizontal  events  below  are  wate' 
crystalline  basement  are  marked  with  arrows  ^op  basement  ref'--'  column  multiples 


OG  -  13 


Fig  4  Details  of  migrated  crossing  profiles  OG-13  (B)  and  OG-8  with  NW  on  profile  OG-8  Some  of  the  basement  faults  coincide  with  offsets  in  the 
processing  details  in  Table  2  Events  interpreted  as  faults  m  the  crystalline  sediment  layering  Fa'ilt  planes  can  be  correlated  between  the  two  n-r,c;9i„n 
basement  are  dipping  mainly  towards  the  SW  on  profile  OG-13  and  towards  profiles 


OG  - 


prominent  basement  faults  the  Mesozoic/Cenozoic  sediments  and  the  seabed 


Fig.  6  The  mapped  faults  on  the  migrated  strike  lines  OG-12,  OG-9,  OG-8 
and  the  crossing  profile  lines  OG-13,  OG-3  and  OG-2  The  line  drawings  are 
time-to-depth  converted  using  an  estimated  velocity-depth  function.  Thick 
lines  mark  the  most  prominent  faults  on  acount  of  continuity  and  deth 
penetration 
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OG  -  13 


Fig.7  Parts  of  modelled  zero-offset  diffraction  hyperbolas  (broken  lines)  inserted  on  a 
OG-13  seismic  section  (left  segment  of  Fig  4).  Velocity-depth  functions  derived 
from  N-Skagerrak  profiling  lines  (Lie  ct  al,  1992b).  The  hypothesis  that  our 
inferred  faults  (.solid  lines)  should  represent  diffraction  hyperbola  artefacts  is  not 
considered  tenable. 
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TWO-WAY-TIME  SECONDS 


An  exercise  in  automating  seismic  record  analysis  and 
network  bulletin  production 

B.O.  Ruud*\  C.D.  Lindholm^^  and  E.S.  Husebye^^ 

Dept,  of  Geology,  Oslo  llniv.,  Oslo,  Norway 
NTNF/NORSAR,  Kjeller,  Norway 


Abstract 

A  long-standing  problem  in  observational  seismology  is  that  of  automating  network  oper¬ 
ation.  In  this  study  we  report  an  experiment  by  which  we  use  the  detector  described  in 
Ruud  and  Husebye  (1992)  for  automatically  picking  P-  and  S-arrivals  in  local  event 
records  stemming  from  the  Norwegian  Seismic  Network.  For  automatic  epitcenter  deter¬ 
mination  a  novel,  robust  grid-search  method  well  suited  for  estimation  problems  with  non- 
Gaussian  observational  errors  is  introduced  in  order  to  handle  outliers.  Even  several  large 
arrival  time  outliers  did  not  prevent  solutions  close  to  the  “true”  epicenter.  In  our  experi¬ 
ment,  38  local  events  from  the  August  1991  bulletin  were  located.  The  number  of  detect¬ 
ing  stations  varied  from  3  to  10  out  of  a  total  of  1 5  stations.  P-  and  S-picking  errors  were 
small,  mostly  within  0.5  sec  for  both  P  and  S.  Phase  identification  errors,  causing  severely 
wrong  P-  and  S-arrivals  were  more  frequent.  Decent  epicenter  determinations  were 
obtained  even  for  events  with  30-50  per  cent  outliers.  The  average  location  difference  of 
our  “automatic”  solutions  compared  to  those  in  the  manually  produced  bulletin  were  about 
15  km. 


1  Introduction 

The  most  outstanding  feat  in  observational  seismology  was  the  creation  of  the  Jef- 
freys-Bullen  tables  in  1937.  This  was  achieved  on  the  basis  of  eminent  analysis  of  analog 


41 


records  from  low-gain,  mechanical  seismometers  in  combination  with  novel  statistical 
schemes  for  event  location  and  travel  time  analysis  (Jeffreys,  1932, 1967).  Such  an 
achievement  would  hardly  be  possible  today,  given  the  same  type  of  recordings,  in  view  of 
the  deteriorating  state  of  routine  seismogram  analysis  at  many  observatories  around  the 
world.  In  an  effort  to  correct  this  situation,  the  seismological  community  through  lASPEI 
has  launched  the  ISOP-project  (International  Seismological  Observation  Period)  under 
NEICAJSGS  management  for  re-establishing  past  skills  in  seismogram  analysis  and 
reportings  to  international  agencies  like  NEIC  and  ISC. 

The  basic  role  of  the  seismological  observatories  is  to  provide  data  for  all  kinds  of 
seismological  research,  and  an  important  product  in  this  respect  is  the  daily  event  bulletin. 
The  mentioned  deterioration  of  routine  seismogram  analysis  is  tied  to  a  curtailing  of  fund¬ 
ing  for  such  work  while  the  workload  is  ever-increasing  due  to  deployment  of  more  and 
more  modem,  high-gain  digital  seismometers.  As  in  similar  situations,  a  viable  solution 
might  be  to  automate,  to  the  extent  possible,  all  aspects  of  the  seismogram  analysis  includ¬ 
ing  the  bulletin  production.  This  is  the  problem  dealt  with  in  this  paper,  and  our  observa¬ 
tional  data  stem  from  the  Norwegian  Seismograph  Network  (NSN)  operated  by  the 
University  of  Bergen.  We  note  that  several  novel  developments  in  seismic  network  opera¬ 
tions  have  been  published  in  recent  years(Bache  et  al,  1990;  Bratt  et  al,  1990;  Boschi  et  al. 
1990,  1991;  Masse  and  Person,  1991). 

2  Network  configuration,  instrumentation  and  data  acquisition 

The  NSN  configuration  is  shown  in  Fig.  1.  TTie  network  Hub  is  physically  located  at 
the  University  of  Bergen  (station  BER)  where  the  routine  data  processing  takes  place. 
Note  that  KONO  and  the  NORESS  and  ARCESS  arrays  do  not  constitute  an  integral  part 
of  the  NSN,  although  their  detection  logs  and  associated  waveform  data  are  accessible.  A 
comprehensive  description  of  the  field  data  acquisition  is  given  elsewhere  (Havskov  et  al. 


42 


1992),  so  here  only  the  essential  features  of  the  system  are  given.  The  NSN  is  divided  in  a 
southern  and  a  northern  segment  along  latitude  62°N  (Table  2).  The  mam  difference  in 
network  operation  is  that  the  southern  stations  transmit  their  recordings  continuously  in 
analog  form  to  the  HUB  via  leased  telephone  lines,  while  for  the  northern  stations  a  simple 
STA/LTA  detector  is  run  on  the  site  and  only  event  recordings  are  transmitted  in  digital 
form  to  the  HUB.  Correct  timing  is  ensured  by  radio  clocks  at  each  independent  station 
(accuracy  1  ms).  The  principal  features  of  the  data  acquisition  are  as  follows  (this  applies 
both  to  field  and  Hub  processes): 

•  The  seismometer,  Kmemetrics  SS-1  Ranger,  “raw”  trace  is  subject  to  an  analog,  5-10 
Hz  bandpass  filter  ( 1 2  dB  per  octave  roll-off)  and  then  transformed  to  digital  form  (50 
Hz  sampling  rate)  via  a  12  bits  A/D  converter  with  an  equivalent  dynamic  range  of 
72  dB. 

•  Signal  screening  via  a  conventional  STA/LTA  detector;  STA  window  is  4  sec;  LTA 
recursively  updated  and  threshold  setting  varies  from  3-6  units.  This  initial  signal 
detection  is  performed  only  on  one  vertical  trace  for  each  station  (irrespective  of 
whether  a  station  has  3C  instrumentation  or  a  cluster  of  sensors).  System  status 
parameters  are  also  exnacted  (Havskov  et  al,  1992). 

•  If  triggering  occurs,  detection  information  (STA/LTA,  trigger  time  and  duration)  plus 
corresponding  waveform  data  are  stored  in  specific  files.  A  ring  buffer  contains  20-45 
hours  (dependent  on  disk  capacity  and  numbers  of  sensors)  of  continuous  recordings 
of  data  for  all  seismometer  traces  at  a  station. 

•  Event  definition  and  subsequent  data  retrieval  is  handled  by  the  Hub  computer. 

The  design  of  the  field  data  acquisition  system  dates  back  to  1985/86.  Planned 
upgrade  will  include  16  bits  A/D  converters  with  gain  ranging  (140  dB),  increased  CPU 


power  and  disk  storage  in  the  field  and  more  stations  will  have  3-component  instrumenta¬ 


tion. 


3  Hub  network  detection  screening  --  event  data  retrieval 

Bihourly  the  Hub  retrieves  detection  and  status  information  from  the  individual  sta¬ 
tions  in  the  network.  Presuming  no  faulty  instrumentation,  the  event  definition  task  is 
resolved  in  the  following  way  (see  Havskov  et  al.  for  more  details); 

•  Minimum  3  stations  triggered  within  a  1 20  sec  time  window  for  either  the  southern  or 
the  northern  NSN  subnets.  For  each  declared  event,  the  retneved  record  segments  start 
30  to  60  sec  prior  to  the  earliest  station  signal  detection  time  and  ends  30  to  60  sec 
after  the  latest  detection  state  cancellation. 

•  The  Hub  retrieves  the  waveform  files  of  uniform  length  from  the  ring  buffers  for  all 
SP  components  for  all  stations  in  the  field  (backup  files  are  used  if  ring  buffer  is  over¬ 
written).  All  waveform  data  is  merged  into  an  event  file  which  constitute  the  base  for 
the  automated  record  analysis  and  the  subsequent  epicenter  determination. 

•  The  system  is  extensively  parameterized,  leading  to  a  high  degree  of  flexibility. 

On  a  daily  basis,  an  average  of  4-6  seismic  events  are  declared,  including  quarry  blasts 
and  mining  explosions. 

4  Automatic  processing  of  event  waveform  data  at  the  Hub 

We  have  here  adapted  the  signal  detector  procedures  detailed  in  Ruud  and  Husebye 
(1992),  hereafter  referred  to  as  Paper  1 ,  for  processing  the  single  vertical  component  and 
the  3C  station  records  of  the  NSN  network.  The  process  is  parameterized,  and  thus  provid¬ 
ing  some  flexibility  in  most  of  the  parameters. 
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•  Signal  quality  check:  Stations  exhibiting  relatively  weak  signals  are  deleted  from  fur¬ 
ther  analysis  according  to  a  simple  criterion.  The  maximum  and  minimum  powers  in  a 
5  sec  moving  time  window  are  computed  and  if  the  max/min  ratio  is  less  than  3  the 
station  is  skipped. 

•  Signal  parameter  extraction:  The  records  are  prefiltered  in  3  passbands,  parameter 
setting  details  in  Table  1 .  For  erch  of  these  an  STA/ITA  detector  is  applied  except  for 
the  first  10  sec  which  is  used  to  initialize  the  LTA  values.  For  3C  stations  P-wave 
polarization  statistics  are  used  in  combination  with  the  STA/LTA  power  detector,  and 
P-slowness  parameters  are  also  estimated  (Roberts  et  al,  1989;  Paper  1). 

•  Signal  parameter  screening  --  phase  identification:  The  detections  within  the  different 
passbands  are  compared  and  the  surplus  (overlapping)  ones  are  deleted  with  basis  in  a 
set  of  decision  rules  (for  details  see  Paper  1).  The  earliest  detection  is  presumed  to  be 
a  P  wave  while  the  second  one,  the  S  wave  (Sg  or  Lg),  is  associated  with  the  strongest 
phase  detection  in  the  time  interval  3  to  100  sec  after  P  onset  (illustrations  in  Figs.  2 
and  3). 

•  Coda  duration  -  local  magnitude  estimation:  The  duration  is  reckoned  from  P  onset 
time  and  until  the  subsequent  coda  power  level  in  a  window  of  length  5  sec  is  less  than 
a  factor  of  1 .6  of  the  initial  LTA  power. 

Note  that  for  local  events  P-S  differential  arrival  times  provide  accurate  epicenter  dis¬ 
tance  estimates,  and  this,  in  combination  with  back-azimuth  estimates  from  individual  3C 
stations,  gives  decent  epicenter  locations  (e.g.,  see  Ruud  et  al,  1988;  Li  and  Thuber,  1991; 
Paper  1 ). 

Automatic  signal  parameter  extraction  is  not  error-proof,  and  with  basis  in  Paper  1 
analysis  we  estimate  that  the  probability  of  picking  both  P  and  S  phases  correctly  at  a  sin- 
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gle  station  to  be  around  0.8.  For  the  NSN-network  (Fig.  1 )  there  are  typically  3-7  reporting 
stations  for  individual  events,  which  implies  that  less  than  30  per  cent  of  the  events  would 
have  “correct”  P  and  S  phase  pickings  at  all  detecting  stations.  Gross  errors  here  could 
adversely  affect  the  event  location  and  the  corresponding  arrival  time  residuals  might  be 
“smeared”  onto  the  other  stations  in  the  network  and  thus  hide  the  erroneous  phase  picks. 
This  so-called  circular  problem  —  a  good  event  location  requires  that  outliers  have  been 
removed  while  outliers  cannot  be  properly  identified  unless  we  have  a  good  location  —  has 
to  be  resolved  to  ensure  automatic  network  operation. 

5  Evaluation  of  the  extracted  signal  parameters 

The  parameters  used  for  epicenter  locations  at  local  distances  are  P  and  S  arrival  times 
and  the  P  slownesses  (that  is,  back-azimuth)  for  individual  stations.  Before  attempting  epi¬ 
center  determinations,  we  undertook  an  appraisal  of  the  quality  of  the  signal  parameters 
extracted  automatically.  In  this  respect  38  events  (Table  3)  from  the  NSN  bulletin  for 
August  1991  were  analyzed.  The  “selection”  criteria  used  were  that  an  epicenter  location 
existed  in  the  manual  bulletin  and  that  digital  records  from  at  least  3  stations  were  at  hand 
for  the  automatic  analysis.  Details  on  network  performance  is  given  in  Tables  2  and  3  and 
in  Fig.  4.  and  5. 

Arrival  time  pickings 

Various  kinds  of  errors  were  found  in  the  automatically  extracted  signal  parameters. 
The  most  frequently  occurring  ones  were  those  due  to  interference  with  “alien”  events  or 
signals  stemming  from  very  local  blasting.  The  arrivals  could  be  picked  correct  in  time  but 
assigned  wrong  phase  identification.  For  very  local  events  it  is  often  difficult  to  separate  P- 
and  S-phases  (A  <  40kni ),  so  only  P  would  be  reported.  In  Figs.  2  and  3  the  above-men¬ 
tioned  problems  are  illustrated  in  terms  of  waveform  displays  including  our  automatically 
picked  P-  and  S-phase  arrival  times.  'ITie  relative  frequency  of  occurrence  of  such  errors 
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can  be  judged  from  Tables  2  and  3.  For  example,  for  Event  1  the  (manual)  NSN  bulletin 
gives  4  P  readings  and  3  S-readings,  while  the  automatic  phase-picker  indicated  5  P-  and  4 
S-readings.  However,  during  the  location  process,  we  find  that  2  i  and  2  S-readings  from 
the  automatic  pickings  were  wrong  (rejected)  as  indicated  in  the  RP  and  RS  columns  to 
the  right  in  Table  3.  It  also  happened  that  the  detector  picked  all  P-  and  S-phases  for  a 
given  event  correctly;  for  our  event  population  this  happened  fori  2  out  of  38  events  (c.  30 
per  cent).  Also,  station  reporting  frequency  varies  considerably  as  indicated  in  Table  2. 
For  example,  the  northern  stations  report  mostly  events  in  the  Kiruna  mining  district  in 
northern  Sweden  (67.7°N,  20.2°E)  (Figs.  3  and  6). 

Automatic  P  and  S  phase  arrival  time  picking  results  are  displayed  in  Fig.  4.  The 
agreements  between  analyst  and  automatic  pickings  are  very  good,  in  particular  for  P 
phases.  There  are  few  outliers,  as  wt  have  manually  corrected  for  phase  identification 
errors  of  the  kinds  discussed  above.  For  S  pha.ses  the  residual  average  is  non-zero,  that  is, 
positively  biased  a  few  tenths  of  a  second.  This  reflects  late  triggering  due  to  increased 
LTA  after  P-onset.  We  consider  the  automatically  picked  P  and  S  arrival  times  to  be  satis¬ 
factorily  precise  compared  to  analyst  pickings,  and  hence  accurate  event  locations  should 
be  attainable. 

Azimuth  estimation 

The  NSN-network  had  5  operative  3C  stations  during  August  1991  (Fig.  1).  For  these 
stations  the  detector  reported  back-azimuths  of  assumed  P-phases  plus  a  phase  identifica¬ 
tion  measure,  the  so-called  PS-paramelei  (details  in  Paper  1 ).  This  parameter  indicates 
whether  the  particle  motion  is  P-type  (PS  =  10  -  pure  P)  or  S-type  (PS  =  0  —  pure  S).  In 
Tables  2  and  3  and  Fig.  5  the  performance  of  azimuth  observation  is  indicated.  The  accu¬ 
racy  of  the  accepted  azimuth  estimates  at  the  NSN  station  were  close  to  3C  azimuth  obser¬ 
vations  from  NORESS,  ARCESS  and  Finnish  stations  (Paper  1;  Suteau-Henson,  1991; 
Tarvainen,  1992)  and  could  apparently  be  improved  by  correcting  for  a  bias  at  some  sta- 
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tions.  However,  many  P-wave  azimuth  determinations  failed  to  meet  our  PS-criterion, 
(Table  1).  Estimation  stability  would  probably  improve  with  instrumentation  having  rela¬ 
tively  better  dynamic  range  and  bandwidth.  Also  scattering  from  the  rough  topogr^hy 
near  the  NSN  stations  (compared  to  for  instance  NORESS  and  ARCESS)  might  corrupt  P- 
wave  polarization  properties. 

Signal  amplitude  (STA)  and  coda  duration  measures 

The  detector  extracted  both  these  parameters  for  each  event  detected.  For  3C  stations 
the  STA  ratio  between  different  components  is  rather  unstable  and  was  not  found  useful 
for  phase  identification.  In  contrast,  the  STA  of  S-phases  is  almost  always  larger  than  STA 
of  P  for  local  events.  The  coda  duration  parameter  as  estimated  was  rather  unstable;  small 
disturbances  in  the  late  coda  could  change  their  value  by  a  factor  of  2.  This  parameter  was 
therefore  not  used  in  magnitude  estimation.  Magnitude  estimation  based  amplitude  (time 
or  frequency  domain)  is  feasible,  but  not  yet  implemented.  However,  we  consider  local 
magnitude  (Ml)  estimation  an  important  task,  which  can  be  resolved  (Alsaker  et  al.,  1990; 
BSth,  1981 ;  Michaelson,  1 990).  For  example,  the  Ml  parameter  may  serve  as  an  important 
measure  for  flagging  local  events  of  marginal  interest,  that  is,  the  numerous  small  (Ml  < 

1 .0)  man-made  chemical  explosions. 

Comments 

The  outcome  of  our  analysis  of  2(X)  station  records  from  the  NSN-network,  displayed 
in  Tables  2  and  3  and  Fig.  4  and  5  are  considered  gratifying.  In  most  cases  4  or  more  reli¬ 
able  signal  parameters  would  be  available  for  epicenter  locations.  This  would  suffice  for 
estimating  the  3  unknowns  origin  time,  latitude  and  longitude  for  a  fixed  focal  depth. 
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6  Robust  earthquake  location  methods  and  results 

For  a  long  time  it  has  been  known  that  errors  in  earthquake  arrival  times  follows  a 
non-Gaussian  distribution  (Jeffreys,  1932).  Consequently,  least  square  estimation  tech¬ 
niques  (e.g.,  Geiger,  1910  and  numerous  modifications)  cannot  yield  optimal  hypocenter 
estimates  (Buland,  1986).  Still,  least  square  algorithms  remain  the  most  commonly  used 
method  for  locating  seismic  events.  The  reason  for  this  is  probably  that  the  theory  behind 
least  squares  is  well  known  and  that  the  method  is  computationally  attractive.  However, 
when  applied  to  non-linear  problems,  convergence  is  not  necessarily  ensured  unless  a 
good  trial  solution  is  known. 

Anderson  (1982)  considers  a  class  of  estimators  particularly  suitable  for  data  sets  with 
large  errors.  He  suggested  an  iterative  algorithm  with  weights  dependent  on  “current” 
residuals  to  minimize  the  objective  function.  Sambridge  and  Kennett  (1986)  and  Nelson 
and  Vidale  (1990)  also  investigated  robust  hypocenter  location  methods  but  used  a  grid- 
search  technique  to  find  the  minimum  of  the  objective  function.  Besides  that  grid-search 
methods  are  better  suited  for  highly  non-linear  optimization  problems,  they  are  also  useful 
for  illustrating  the  “resolution”  of  the  estimated  epicenter.  Buland  (1976)  made  contour 
maps  of  the  sum  of  squared  residuals  and  also  derived  their  relation  to  corresponding  con¬ 
fidence  regions  assuming  Gaussian  errors.  Prugger  and  Gendzwill  (1988)  produced  simi¬ 
lar  maps  using  average  absolute  residuals.  They  recommended  using  the  Simplex  step 
algorithm  (Nelder  and  Mead,  1965)  to  solve  the  non-linear  minimization  problem. 

In  our  experiment,  back-azimuth  estimates  were  available  from  some  3-component 
stations  so  we  also  wanted  to  include  these  in  the  location  scheme,  that  is,  in  addition  to 
the  P-  and  S-readings.  The  input  data  for  our  location  algorithm  thus  consist  of  different 
kinds  of  physical  quantities  so  it  was  necessary  to  define  normalized  residuals  by  dividing 
by  a  priori  estimates  of  the  standard  error  of  each  observation,  i  .e.. 
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do-  -  dc- 

where  dOj  and  dcj  are  observed  and  calculated  data  respectively,  and  a.  is  the  correspond¬ 
ing  standard  error.  The  calculated  residuals  are  functions  of  the  hypocenter  coordinates 
(latitude,  longitude,  depth  and  origin  time),  the  coordinates  of  the  observing  station  and 
the  kind  of  observation  (azimuth,  or  in  case  of  arrival  times  a  phase  designation,  see  exam¬ 
ples  in  Table  4). 

Initial  experiments  with  both  L2  (least  squares)  and  L|  (absolute  errors)  norms 
showed  that  a  few  large  errors  could  heavily  bias  the  epicenter  estimate,  i.e.,  pull  it  away 
from  the  “true”  epicenter.  Our  next  step  was  to  try  an  objective  function  of  the  form 


/(r) 


(as  compared  to  for  the  L2  norm  and  norm).  Here  Wj  axe  weights 

that  depend  on  the  confidence  we  have  in  each  observation  and  the  summation  is  over  all 
observations  (both  arrival  times  and  azimuths  for  all  stations).  This  estimator  is  very  simi¬ 
lar,  although  not  identical,  to  some  of  the  estimators  considered  by  Anderson  (1982).  The 
advantage  of  this  objective  function  is  that  it  is  practically  fiat  (insensitive)  for  |r.|  >  >  1 , 
while  for  <  <  1  it  behaves  like  the  1-2  norm. 

The  hypocenter  location  problem  now  becomes  very  non-linear  as  compared  to  least 
square  methods  and  a  grid  search  approach  is  necessary  to  find  the  global  minimum  of  the 
objective  function.  As  most  of  the  events  were  explosions  and  the  remaining  earthquakes 
were  relatively  shallow  (within  the  crust),  we  used  a  fixed  zero  focal  depth  and  searched 
over  a  gnd  of  latitude  and  longitude  lines.  An  initial  set  of  5 1  x  51  points  covering  an  area 
of  about  1000  X  KXK)  km  surrounding  the  detecting  stations  proved  sufficient  for  all  the 
local  events.  For  each  grid  point  an  additional  search  over  time  was  performed  to  find  the 
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best  origin  time.  Typical  execution  times  were  about  I  minute  per  event  on  a  SUN  Sparc 
station.  A  simple  crustal  model  consisting  of  two  horizontal  layers  over  a  halfspace  was 
used  with  travel  times  calculations  for  every  point  in  the  grid.  Tabulated  travel  times  could 
be  used  to  save  time,  especially  if  more  complex  models  were  to  be  considered  (see  Nel¬ 
son  and  Vidale,  1990).  In  the  location  experiment  we  presumed  that  the  automatically 
picked  P-  and  S-phases  were  the  first  arriving  phase  of  each  kind.  Specific  secondary 
phases  like  Pg  and  Lg  (at  distances  beyond  200  km)  are  sometimes  observed  and  can  be 
allowed  for  in  the  location  algorithm,  but  were  not  used  in  the  automatic  processing.  Our 
observational  data  indicated  that  P-pickings  were  more  reliable  than  S-pickings,  so  that 
the  former  were  given  greater  weight  (2.0  for  P;  1.0  for  S  and  P-azimuth).  The  a  priori 
standard  deviations  were  set  rather  high  to  compensate  for  the  coarse  initial  grid:  2.0  sec 
for  P,  3.0  sec  for  S  and  15  deg  for  the  P-azimuth. 

The  results  presented  here  (Table  3)  were  computed  with  the  rather  coarse  initial  grid, 
but  if  desiered  a  more  refined  “final”  hypocenter  estimate  could  be  found  in  several  ways: 
(i)  identify  outliers  given  a  cutoff  criteria  and  proceed  with  one  of  the  commonly  used 
location  programs  (HYP071  or  variants  hereof,  e..g.,  Lienert  et  al,  1986),  (ii)  restart  the 
search  using  a  denser  grid,  or  (lii)  proceed  with  an  iterative  algorithm  and  the  same  objec¬ 
tive  function  (convergence  should  not  be  a  problem  when  sufficiently  close  to  the  global 
minimum). 

Practical  demonstrations  of  the  above  epicenter  location  method  are  illustrated  in 
Figs.  6  and  7  and  Table  4  for  the  two  events  in  Figs.  2  and  3.  The  full  scale  test  here  was  to 
locate  38  events  using  the  signal  parameters  automatically  derived  at  the  Hub.  The  differ¬ 
ences  between  the  analyst  derived  locations  (H YP07 1 )  and  ours  have  a  standard  deviation 
of  18  km  in  the  south-north  direction  and  14  km  in  east- west  (Table  3). 
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7  Discussion 


We  have  here  demonstrated  automatic  network  operation  in  the  sense  that  a  daily  bul¬ 
letin  can  be  produced  which  would  not  be  much  inferior  to  that  of  an  analyst.  With  the 
data  used  the  success  rate  of  our  automatic  event  analyser  to  extract  reliable  P-  and  S- 
arrivals  was  around  70-80  per  cent  for  each  station.  This  entails  that  only  around  30  per 
cent  of  the  processed  events  would  have  correct  arrival  times  and  phase  identification  for 
all  contributing  stations. 

With  more  operational  experience  (and  possibly  system  upgrading),  the  performance 
of  the  automatic  analysis  is  likely  to  be  improved.  However,  as  we  cannot  safeguard 
against  gross  observational  errors,  a  robust  epicenter  location  scheme  is  required.  The 
method  introduced  here  fulfills  such  requirements  even  for  rather  adverse  cases  where 
faulty  observations  could  reach  nearly  50  per  cent  of  the  total  (see  Table  2). 

The  main  results  from  this  study  is  condensed  in  Figs.  5-7  and  Table  2  and  3,  but  we 
would  like  to  give  a  few  comments  on  the  somewhat  problematic  mining  explosions  in  the 
Kiruna  area.  Firstly,  Event  19  is  interesting  in  view  of  the  relatively  large  mismatch 
between  the  bulletin  solution  in  the  northern  Bothnian  Bay  and  ours.  A  comparison  with 
the  Helsinki  bulletin  gave  that  the  latter  was  approximately  correct,  that  is,  a  Kiruna  min¬ 
ing  explosion.  The  bulletin  mislocation  was  simply  due  to  the  use  of  the  Pg  notation  for 
some  of  the  first  arrivals.  Another  somewhat  problematic  Kiruna  event  was  no.  23,  where 
neither  the  bulletin  nor  our  solution  is  particularly  good.  The  advantage  of  the  latter  is  that 
with  our  scheme  for  epicenter  determinations  we  can  obtain  a  visual  measure  of  the  epi¬ 
center  resolution,  as  demonstrated  in  Figs.  6  and  7.  In  the  latter  case  the  resolution 
becomes  relatively  poor,  particularly  in  the  N  W-SE  direction  (Fig.  6).  In  general,  the  NSN 
network  geometry  is  far  from  ideal  because  it  is  (for  geographical  reasons)  restricted  to  a 
very  elongated  area.  For  accurate  locations  of  the  many  events  outside  the  network  area  it 
IS  crucial  to  have  good  S-phase  readings  (Buland,  197(>;  Comberg  el  al,  1990). 
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The  events  in  Table  3  are  local  m  the  sense  that  all  station-epicenter  distances  are  less 
than  750  km.  In  cases  of  regional  and  teleseismic  events  the  above  location  scheme  would 
be  inconvenient,  and  should  instead  be  used  in  combination  with  those  presented  by,  for 
example,  Bratt  and  Bache  ( 1 988),  Cassidy  et  al  ( 1 989)  and  Ruud  ( 1 990). 

Finally,  we  want  to  discuss  our  results  here  m  the  context  of  generalized  seismological 
network  operation.  Firstly,  in  the  current  NSN  operation  event  detection  is  a  two-step  pro¬ 
cess.  A  rather  crude  STA/LTA  detector  is  used  in  the  field,  while  the  detector  described  in 
Paper  1  is  used  for  extracting  signal  parameters  from  segmented  waveform  data  at  the 
Hub.  The  latter  detector  could  naturally  be  installed  at  the  individual  stations  in  the  field 
given  sufficient  CPU  capacity.  The  advantage  here  is  that  for  very  many  events  it  would 
suffice  to  transmit  on  a  routine  basis  only  parameterized  (signal  parameters)  event  data  to 
the  Hub. 

We  have  earlier  discussed  a  few  “restart”  schemes  for  more  refined  event  locations. 
Better  results  here  would  be  obtained  if  such  undertakings  were  combined  with  reprocess¬ 
ing  and  reana'ysis  of  the  original  'ecords  (e.g.,  see  Roberts  and  Christoffersson.  1 9*^0).  For 
example,  given  ,  .liminary  epicenter  coordinates  this  information  would  facilitate  pick¬ 
ings  of  individual  P-  and  S-phases  like  Pg,  Pn.  Sg  (Lg)  and  Sn,  which  in  turn  could  help  us 
to  provide  decent  focal  depth  estimates  (Ruud  et  al,  1988;  Thurber,  1989).  Further  refine¬ 
ments  would  be  feasible  if  access  to  a  priori  seismicity  information,  similar  to  the  knowl¬ 
edge-based  analysis  system  introduced  for  arrays  (Bache  et  al,  1990;  Bratt  et  al,  1990).  An 
illustrative  example  here  is  the  time-space  stability  of  many  of  the  Kiruna  mining  explo¬ 
sions  (events  3,8,15,23,27,31  and  37  in  Table  3). 

The  NSN  network  comprises  relatively  few  stations  (Fig.  1),  while  in  operational 
terms  we  distinguish  between  a  southern  and  a  northern  subnet.  The  rationale  here  is  two¬ 
fold:  namely,  i)  to  facilitate  the  event  association  problem  and  ii)  little  overlap  in  event 
populations  between  the  two  subsets.  In  the  latter  case  another  example  is  that  the  event 
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overlap  between  the  NORESS  and  ARCESS  arrays  (Fig.  1)  is  only  about  10  per  cent.  The 
reason  for  this  is  that  the  numerous  man-generated  events  are  of  ML  magnitudes  less  than 
2.0,  and  thus  in  general  are  not  recorded  at  epicenter  distances  beyond  500  km  (Paper  1 ; 
Suteau- Henson,  1991).  Regarding  event  association,  that  is,  the  problem  of  determining 
which  of  the  delected  phases  stem  from  a  specific  event,  its  severity  increases  with 
increasing  station  separations.  For  network  subsets  of  aperture  3-500  km,  event  associa¬ 
tion  is  not  much  of  a  problem.  For  netwoiks  much  larger  than  the  NSN  subdivision  into 
several  subsets  of  the  mentioned  size  could  well  be  done. 

There  is  a  growing  interest  to  utilize  3C  data  more  extensively  in  the  seismological 
community.  It  is  therefore  interesting  to  investigate  the  importance  of  the  additional  infor¬ 
mation  provided  by  the  3C  stations  compared  to  the  single,  vertical  component  (1C)  sta¬ 
tions.  In  this  study  we  could  not  fully  lest  performance  differences  between  3C  and  IC 
stations  since  the  analysed  data  were  restricted  to  the  NSN  event  selection  based  on  1C 
station  detection.  For  this  experiment  most  P-  and  S-arrival  times  could  be  picked  rela¬ 
tively  accurate  for  both  type  of  stations.  In  a  “true”  detector  environmeml,  with  many  low 
SNR  events,  we  would  expect  3C  stations  to  have  better  detection  and  time  picking  abili¬ 
ties  than  a  single  component  station  due  to  lower  threshold  for  P-phase  detection  and 
because  S-phase  often  are  stronger  on  the  horizontal  components. 

Regarding  the  use  of  azimuth  information  we  found  that  when  a  sufficiently  large 
number  of  stations  and  correct  arrival  limes  (at  least  3  stations  and  5  arrival  times)  were 
available  the  azimuth  information  did  not  contribute  much  to  the  epicenter  location.  How¬ 
ever,  azimuth  is  a  crucial  parameter  when  locating  events  for  a  single  station  and  is  also 
very  important  with  two  detecting  stations.  For  high-frequency  “local”  signals,  P-phases 
are  not  necessarily  well  polaiized  and  in  many  cases  the  3(!  stations  failed  to  provide  azi¬ 
muth  estimated  (Table  2).  Also  with  respect  to  P/S  phase  drscrimination  the  results  from 
the  NSN  network  were  less  impressive  than  earlier  results  at  NORliSS  (Paper  1).  The 
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most  likely  reason  for  the  deteriorated  P-wave  polai  ization  is  scattering  (from  P  to  surface 
waves)  due  to  the  very  rugged  topography  near  most  of  the  NSN  stations.  This  illustrates 
that  polarization  information  can  be  difficult  to  utilize  m  an  on-line  detection  environment. 

The  full  advantage  of  3C  recordings  is  probably  better  realized  in  a  post-processing 
stage  when  a  preliminary  event  location  is  available.  When  an  approximate  arrival  time  is 
known  more  complex  particle  motion  models  can  be  introduced  for  wavefield  decomposi¬ 
tion  analysis.  This  could  enable  us  to  identify  secondary  arrivals  which  are  important  for 
focal  depth  estimation. 


8  Concluding  remarks 

The  results  of  this  study  implies  that  seismic  network  operation  in  terms  of  daily  bul¬ 
letin  generation  to  a  large  extent  can  be  automated.  This  would  enable  the  analyst  to  con¬ 
centrate  on  more  in-depth  analysis  of  real  earthquake  recordings,  which  is  a  principal  goal 
of  the  ISOP  project.  Today,  it  appears  that  far  too  much  analyst  effort  is  spent  on  analysis 
of  local  explosions  which  are  of  marginal  scientific  value. 

Strictly  speaking,  the  conclusions  drawn  here  are  restricted  in  the  sense  that  they  are 
based  on  our  experience  with  recordings  from  the  Norwegian  Seismograph  Network 
(NSN).  To  what  extent  our  results  applies  toother  seismic  networks  and  stations  sited  in 
different  tectonic  envnonments  remains  partly  an  open  question  However,  based  on  our 
experience  with  analysis  of  various  data  we  anticipate  that  similar  results  can  be  obtained 
at  most  networks  by  adjusting  the  detector  parameter  settings  and  by  modifications  of  the 
associated  decision  rules. 
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Filter  band  (Hz) 


1-4 

4-8 

8-16 

STA  window  (sec) 

1.00 

0.60 

0.40 

Step  length  (sec) : 

0.25 

0.15 

0.10 

LTA  delay  (sec) 

15.00 

9.00 

6.00 

Threshold_l  : 

1.70 

1.80 

1.90 

Threshold._2 

2.55 

2.60 

2.75 

Pc  min 

0.40 

0.40 

0.40 

K-pararoeter 

4 

4 

4 

Table  1 .  Most  entries  are  self-explanatory  and  otherwise  detailed  in  Paper  1 .  Step  length  or 
STA  updating  interval  is  1/4  of  STA  window  length.  LTA  delay  is  the  time  gap  between 
the  STA  and  LTA  windows.  The  two  thresholds  are  for  linearly  polarized  (1)  and  for  unpo- 
larized  (2)  signals.  Pc_min  is  a  threshold  for  polarization  (predicted  coherence  limit).  K- 
parameter  gives  the  minimum  number  of  consecutive  windows  that  must  be  triggered 
before  the  detection  is  accepted  (duration  limit). 
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- 
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- 

ASKl 

2.*) 
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4 

16 

4 
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1 

21 

1 

- 

- 

HYA 

24 
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18 

2 

- 

- 

FOO 
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0 

6 

0 

6 

1 

FRO 

6 

0 

6 

0 

- 

- 

MOL 

16 

5 

11 

0 

- 

- 

NSS 

4 

1 

3 

2 

- 

- 

MORI 

6 

2 

2 

0 

5 

2 

LOF 

11 

3 

10 

3 

5 

0 

TRO 

10 

0 

9 

3 

8 

0 

KTKI 

10 

0 

10 

3 

- 

•  •  *  r  __  _ 

Table  2.  The  Norwegian  Seismograph  Network  (NSN)  station  listing.  In  the  table  the  col¬ 


umns  OP,  RP,  OS  and  RS  give  numbers  of  observed  (O)  and  rejected  (R)  (during  epicenter 
locaoon)  P-  and  S-pickmgs.  The  same  applies  to  the  two  rightmost  OAzi  and  RAzi  col¬ 


umns  conserning  azimuth  observations. 
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Table  3.  Event  tabulations  and  analysis  results.  The  columns  to  the  left  give  origin  time, 
epicenter  coordinates  and  focal  depths  as  taken  from  the  Bergen  seismological  bulletin  for 
August  1991 .  The  next  three  columns  give  the  number  of  observing  stations  and  the  num¬ 
ber  of  reported  P-  and  S-  arrival  times.  The  slashes  (B/A)  separate  reportings  in  the  bulle¬ 
tin  (B)  and  in  the  automatic  detector  log  (A).  The  NA  column  indicates  the  number  of 
reported  azimuths.  The  next  3  columns  give  the  differences  between  the  bulletin  listed  on- 
gin  times  (dOT)  and  epicenter  coordinates  (dLat;  dLong).  The  3  rightmost  columns  (RP, 
RS  and  RA)  give  the  number  of  P-  and  S-phase  pickings  and  azimuth  estimates  which 
were  rejected  (neutralized)  by  our  epicenter  determination  scheme.  Note,  when  averaging 
the  location  errors,  events  no.  19  and  23  were  deleted  as  obviously  the  bulletin  solutions 
here  are  not  very  good  (text).  Likewise.  dO  L  averaging  was  exclusive  of  events  for  which 
the  bulletin  gave  focal  depth  in  excess  of  3  km  -  we  u.se  zero  depth  as  default  value.  In  the 
bottom  line,  sums  and  averages  are  given  for  various  parameters.  The  standard  deviation 
in  dl^t  and  dLong  were  0.18  and  0.14  respectively. 
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Sum:  NP(A)vNS(A);s'A,RP,RS,RA; 

Average:  dOT, dLal, dLong:  202  145  40  M  0.02  -0.03 


Event  1 9 


Event  2  5 


Origin 

tine  (sec) : 

-19.99 

Latitude 

(deg  N) 

67.10 

Longitude 

(deg  E) : 

20.80 

Fixed 

depth  (kn) 

0.00 

Obs 

Res 

Wei 

NSS 

P 

49.40 

0.55 

2.00 

NSS 

S 

121.90 

22.78 

1.00 

KTKl 

P 

16.50 

-0.30 

2.00 

KTKl 

S 

45.00 

1.35 

1.00 

TRO 

P 

24.00 

0.17 

2.00 

TRO 

S 

58.50 

2.69 

1.00 

TRO 

Az 

170.00 

6.11 

1.00 

LOF 

P 

27.30 

-1.00 

2.00 

LOF 

S 

74.90 

11.35 

1.00 

LOF 

Az 

110.00 

2.93 

1.00 

MOR 

P 

22.40 

-0.24 

2.00 

MOR 

S 

60.10 

6.34 

1.00 

MOR 

Az 

83.00 

15.61 

1.00 

Origin 

tine  (sec) 

43.36 

Latitude 

(deg  N) 

61.05 

Longitude 

(deg  E) 

4.40 

Fixed 

depth  (kn) 

0.00 

Obs 

Res 

Wei 

ODDI 

P 

92.49 

20.65 

2.00 

MOL 

P 

106.70 

26.43 

2.00 

SUE 

P 

47.40 

0.84 

2.00 

HYA 

P 

58.30 

-1.02 

2.00 

HYA 

S 

70.80 

-0.15 

1.00 

BER 

P 

68.78 

10.72 

2.00 

EGD 

P 

28.94 

-32 . 00 

2.00 

EGD 

S 

71.60 

-2.16 

1.00 

ASK 

P 

55.82 

-0.13 

2.00 

ASK 

S 

65.58 

0.45 

1.00 

FRO 

P 

57.30 

0.32 

2.00 

FRO 

S 

67.30 

0.39 

1.00 

FOO 

P 

55.80 

0.92 

2.00 

FOO 

S 

63.60 

0.32 

1.00 

FOO 

Az 

241.00 

31.30 

1.00 

Table  4.  Listing  of  results  from  the  of  automatic  processing  of  events  19  and  25  in  Table  3. 
All  times  are  in  seconds  relative  to  15:39:00  for  event  19  and  17:05:00  for  event  25.  Azi¬ 
muths  are  in  decimal  degrees.  The  column  labels  “Obs”  means  observed  values,  while 
“Res”  gives  the  residuals  (observed  -  calculated).  “Wei”  are  the  weights  assigned  to  each 
observation.  For  event  25  the  outliers  are  ODDI  (P),  MOL  (P),  BER  (P),  EGD  (P),  and 
FOO  (Az),  while  for  event  19  the  NSS  (S),  MOR(S)  and  LOF(S)  are  outliers.  The  same 
events  are  shown  in  Figs.  3  and  4. 
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0°  10'  20'  30° 


0'  10'  20'  30' 

Fig.  1 .  TTie  Norwegian  Seismograph  Network  (NSN)  as  of  August  1991 .  Differences  in 
instrumentation  are  indicated:  Open  triangle  -  NSN  single  or  cluster  (MOR  and  KTK)  of 
vertical  component  stations,  filled  triangle  -  NSN  3-component  statio,  filled  diamond  - 
KONO  broad  band  3-component  instrument,  open  circle  -  NORSAR  large  aperture  array 
of  short  period  mstrumcnis.  filled  circles  -  small  aperture  arrays  (NORESS  and  ARCESS) 
of  short  period  vertical  and  3-componcnt  instruments.  KONO  and  the  NORSAR,  NOR¬ 
ESS  and  ARCESS  arrays  are  not  a  part  of  NSN,  although  their  data  are  accessible. 

f)4 


Fig.  2.  Event  records  from  the  southern  segment  of  the  NSN  --  Event  25  in  Table  3.  The 
time  axis  shows  minutes  and  seconds  after  17:00.  Four  of  the  P-phase  pickings  of  the 
detector  were  rejected  during  the  epicenter  determination,  namely,  ODD,  EGD,  BER  and 
MOL.  This  event  illustrates  typical  problems  in  automatic  phase  pickings.  For  example, 
for  ODD,  BER  and  MOL  the  P-phase  is  missed  by  the  detector  and  as  a  consequence  the 
secondary,  clear  S-amvals  are  denoted  P-phases.  For  EGD,  a  weak  preceding  signal  (out¬ 
side  plot  at  17:05;28.9)  triggers  a  “wrong”  P-amval  while  the  S-phase  is  correctly 
detected  and  identified.  For  some  reason  the  detector  missed  the  S-phase  for  station  SUE. 
Although  four  P-phase  arrivals  were  rejected  from  the  subsequent  epicenter  analysis,  the 
difference  relative  to  the  bulletin  solution  is  small  indeed  (see  Table  3). 
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Fig.  3.  Kiruna  mining  blast  as  recorded  by  stations  in  the  northern  segment  of  the  NSN  — 
Event  19  in  Table  3.  The  time  axis  shows  minutes  and  seconds  after  15:00.  The  P-phase 
pickings  coincide  with  those  of  the  analyst,  while  some  S-phase  pickings  are  too  late  (NSS 
outside  plot  at  15:41 :01 .9).  The  bulletin  epicenter  solution  in  Table  3  is  erroneous  relative 
to  ourL-  and  reflects  the  analyst  decision  to  give  a  Pg-notation  for  all  first  arrivals  and  to 
reject  all  S-phases. 
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P  phase  picking  errors 
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Fig.  4.  Histogram  plots  of  P-  and  S-residuals.  The  time  residuals  are  the  differences 
between  the  detector  picked  phase  onset  times  and  those  picked  by  the  analyst  (in  the  bul¬ 
letin).  Errors  due  to  wrong  phase  identification  have  been  removed.  The  P-residuals  con¬ 
sists  of  173  observations  including  one  outside  the  displayed  interval.  For  S-pickings  there 
are  120  observations  with  10  outliers  (all  with  positive  residuals).  The  P-residuals  are 
small  and  besides  of  the  same  order  as  reported  by  Baer  and  Kradolfer  (1987)  for  similar 
experiments  using  recordings  from  the  Swiss  network,  S-residuals  are  larger  and  exhibit  a 
small  positive  bias  of  around  0.4  sec. 
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Azimuth  Residual  (Degrees) 

Fig.  5.  Azimuth  residuals  for  all  3C  sialions.  The  numbers  to  the  light  is  the  total  number 
of  observations  and  the  number  in  the  midle  is  number  of  observations  with  residual  less 
than  30  degree.  The  residuals  were  calculated  as  the  difference  between  the  observed  azi¬ 
muth  and  the  azimuth  computed  from  the  bulletin  event  locations.  (As  the  residuals  were 
calculated  from  a  larger  set  of  event,  the  total  number  of  observations  is  different  from 
Tables  2  and  3.) 
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Fig.  6.  Epicenter  solution  for  Event  1 9  in  Table  3  and  with  station  records  displayed  in  Fig 
4.  The  geometry  of  the  recording  stations  gives  relatively  poor  resolution  in  the  NW/SE 
direction.  Grid  steps  were  0.03  degrees  in  latitude  and  0.10  degrees  in  longitude.  The 
greyshades  shows  the  value  of  the  objective  function  (eq.  2)  for  fixed  zero  depth  and  for 
the  best  origin  time  at  each  grid  point. 
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Fig.  7.  Epicenter  solution  for  Event  25  in  Table  3  and  with  station  records  displayed  in  Fig 
3.  The  network  geometry  is  relatively  symmetric  for  this  event,  hence  a  more  circular  res 
olution  area.  Grid  steps  were  0.05  degrees  m  latitude  and  0.10  degrees  m  longitude.  The 
greyshades  shows  the  value  of  (lie  objective  function  (eq.  2)  for  fixed  zero  depth  and  for 
the  best  origin  time  at  each  gi  id  point. 
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Abstract 

Using  fundamental  model  Rayleigh  (Rg)  recordings  from  seven  arrays  on  four  conti¬ 
nents,  we  have  explored  structural  properties  in  the  respective  siting  areas  through  inver¬ 
sion  of  Rg  phase  velocity  dispersion  characteristics.  The  models  used  were  one  and  two 
layers  over  a  halfspace  with  shear  velocities  and  layer  thickness  as  unknowns.  In  the  latter 
case,  layer  thicknesses  were  fixed  at  0,5  and  1 .0  km,  respectively.  The  estimated  halfspace 
S-velocities  were  remarkably  conststent  between  the  arrays  with  an  average  value  of  3.55 
±  0.03  kms"' .  Likewise,  in  the  case  of  the  one-layer  model,  the  average  S-velocity  was 
2.87  ±  0.07  kms  '  and  thus  also  rather  consistent.  However,  estimated  layer  thicknesses 
ranged  from  0.12  km  (Yellowknife)  to  1 .6  km  (Alice  Springs).  A  comparison  of  the  one- 
layer  and  two-layer  model  results  implies  that  layer  thicknesses  and  corresponding  shear 
velocity  are  coupled.  For  a  belter  physical  insight  into  Rg  propagation  in  complex  media, 
2D  finite  difference  synthetics  were  computed.  Inhomogeneous  media  were  described  by 
0-th  order  von  Kaimaii  functions  (scIf-Mmilai).  In  the  lattei  case,  S-seailering  wavelets 
would  intefere  with  the  Rg  wavetrain  and  thus  reduce  the  accuracy  of  the  phase  velocity 
measurements,  especially  at  short  periods.  Finally,  excitation  of  Rg  waves  drops  sharply 
for  focal  depths  below  2-3  km. 
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1  Introduction 

In  a  previous  paper  we  discussed  and  demonstrated  the  usefulness  of  using  short- 
period  fundamental  mode  Rayleigh  (Rg)  waves  as  a  tool  for  mapping  upper  crustal  struc¬ 
tures.  'I'he  observational  data  were  local  explosion  recordings  from  15  NORSAR  subai- 
rays  (aperture  ca  7  km)  located  in  SE  Norway  (Lxikshtanov  et  al,  1991).  In  this  study  we 
report  on  a  similar  undertaking  using  local  seismogram  records  from  6  other  arrays 
located  on  4  continents  (Table  1).  Traditionally,  seismic  reflection  and  refraction  profiling 
are  used  for  crustal  mapping,  but  none  of  these  techniques  are  particularly  effective  as 
regards  the  uppermost  2-3  km  of  the  crystalline  crust.  Not  surprisingly,  some  investigators 
have  used  the  Rg  part  of  the  records  for  detailed  mapping  of  the  upper  crustal  layering 
along  such  profiles  (e  g..  MacBcth  and  Burton,  1985,  1986;  Reiter  el  al,  1988;  and  Saikia 
et  al,  1990).  A  drawback  with  this  kind  of  surveys  is  the  high  cost  of  the  field  work  m  con¬ 
trast  to  high-quality  recordings  from  stationary  arrays  which  are  easily  available  and  free 
of  charge. 

In  this  study,  upper  crust  structural  information  bearing  on  the  respective  array  sites 
are  derived  by  inversion  of  frequency-dependent  Rg  phase  velocity  observations.  Compu¬ 
tational  details  on  the  dispersion  analysis  and  inversion  formalism  used  are  given  in  the 
Lokshtanov  et  al  ( 1991 )  paper  and  will  not  be  repeated  here.  Furthermore,  as  we  are  deal¬ 
ing  with  Rayleigh  wave  observations  in  different  geological  environments,  relative  Rg 
excitation  and  propagation  efficiency  are  explored  through  synthetic  seismogram  analysis. 
In  this  respect,  also  heterogeneous  or  random  scattering  media  would  be  considered. 

2  Observational  data,  their  analysis  and  results 

Our  observational  data  stem  from  six  small-to-medium  sized  arrays  on  4  continents, 
namely,  ARCESS  (ARAO.  Norway),  Eskdalemuir  (EKA,  Scotland),  GERESS  (GEC2, 
Germany).  Garibidinaur  (GB  A,  India),  Alice  Springs  (ASAR.  Australia)  and  Yellowknife 
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(YKA,  Canada).  For  the  sake  of  completeness,  we  have  included  the  results  from  the 
NORESS  (NRAO,  S.  Norway)  already  published  by  Ixikshtanov  et  al  (1991).  The  neces¬ 
sary  technical  details  are  given  in  Table  1 ,  including  the  number  of  earthquake/explosion 
recordings  used  for  the  individual  arrays.  Representative  waveform  records  are  shown  in 
F'lg.  1. 

As  is  well  known,  crust  and  upper  mantle  structural  information  can  be  denved  from 
Rayleigh  wave  dispersion  analysis.  The  standard  procedure  here  is  first  to  extract  disper¬ 
sion  information,  that  is,  phase  and/or  group  velocities  and  tl.  n  sol  ve  the  inverse  problem 
in  terms  of  medium  structural  parameters.  The  structure  is  represented  through  a  stack  of 
uniform  layers  and  the  model  parameters  introduced  in  the  inversion  formalism  are  com- 
pressional  and  shear  wave  velocities,  layer  densities  and  thicknesses.  Further  details  in 
Lokshtanov  et  al  (1991),  hereafter  referenced  as  Paper  1 . 

A  conventional  f-k  technique  was  the  principal  tool  in  the  dispersion  analysis.  Signal 
power  was  extracted  as  a  function  of  phase  velocity  for  a  given  frequency  and  azimuth,  as 
illustrated  in  Fig.  2.  The  scatter  in  the  observations  is  quite  moderate,  as  also  experienced 
in  our  Rg  dispersion  analysis  for  the  NORSAR  subarrays  (Paper  1).  Hence,  2-3  events  for 
each  array  appear  adequate  for  dispersion  estimates  which  are  shown  in  Fig.  3. 

Our  inversion  algorithm  is  limited  to  models  with  constant  physical  parameters  within 
horizontal  layers,  while  in  the  real  crust  this  assumption 's  likely  to  be  only  approximately 
valid.  Also,  the  limited  bandwidth  of  the  observations  (typically  0.6  -  1 .8  sec),  the  smooth¬ 
ness  of  the  dispersion  curves  and  the  absence  of  higher  mode  Rg  waves  imply  that  we  can¬ 
not  resolve  more  than  2-3  model  parameters.  Hence  we  have  restricted  our  inversion 
experiments  to  two  models,  each  having  three  free  parameters.  These  unknowns  are  shear 
velocities  and/or  layer  thickness,  while  densities  are  fixed  and  P  and  S  velocities  are  cou¬ 
pled  through  a  Poisson  ratio  of  0.25.  Model  1  consists  of  a  single  layer  over  a  halfspace 
with  the  two  shear  velocities  and  layer  thickness  as  unknowns.  Model  2  consists  of  two 
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layers  (thicknesses  0.5  and  1 .0  km)  over  a  halfspace  and  the  three  shear  velocities  are  the 
free  model  parameters. 

Model  1:  One  layer  over  half  space 

The  ^2  estimates  of  S-wave  velocities  in  the  upper  ciust  (half-space  equivalent)  are 
remarkably  consistent,  varying  from  3.44  kms  '  (YKA)  to  3.68  kms"^  (ASAR)  or  by  7  per 
cent  (details  in  Table  2  and  Fig.  4a).  The  exception  here  is  EKA  at  3.10  kms'^  which  was 
expected  in  view  of  its  siting  on  sedimentary  rocks.  The  Pj  estimates,  representative  for 
the  top  crustal  layer,  are  also  rather  consistent,  varying  from  2.59  kms'^  (EKA)  to  2.96 
kms'*  (ASAR),  or  by  1 3,5  per  cent.  The  top  layer  thicknesses  range  from  0.12  km  (YKA) 
to  1 .6  km  (ASAR)  and  thus  vary  considerably.  The  latter  (ASAR)  is  not  directly  compara¬ 
ble  to  the  other  arrays,  since  the  observations  are  confined  to  relauvely  long  wavelengths. 
Layer  thickness  in  the  range  0.5-1 .0  km  are  most  common  for  the  siting  areas  and  fre¬ 
quency  ranges  considered. 

Model!:  Two  layers  over  lialfspace 

For  this  model  the  layer  thicknesses  were  fixed  at  respectively  0.5  and  1 .0  km,  so  here 
the  unknowns  were  restricted  to  layer  and  half-space  shear  velocities  and  subjected  to  the 
constraint  that  Pj  ^  P2  ^  P,  The  half-space  velocities  (P^  estimates)  are  very  similar,  as 
seen  from  Table2  and  Fig.  4b.  For  those  arrays  with  relatively  thin  Model  1  layer  thick¬ 
nesses  (ARAO,  EKA,  GBA  and  YKA),  the  layer  2  velocities  (P2  estimates)  are  less  than 
10  per  cent  below  their  respective  half-space  velocities.  For  NR  AO  and  GBA  having  rela¬ 
tively  large  Model  1  layer  thicknesses,  the  Pj  velocities  have  been  lowered  in  comparison 
to  Model  1 .  The  overall  model  fit  parameter  a  (also  given  in  Table  1 )  is  small  for  all 
arrays  for  both  Model  1  and  2,  thus  implying  a  good  fit  between  observations  and  derived 
model  parameters.  However,  model  resolution  is  somewhat  poor  in  particular  for  ASAR, 
which  in  turn  reflects  lack  of  short  wavelength  observations. 
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3  Synthetic  Rg  modelling 

The  Models  1  and  2  above  are  homogeneous  models  in  the  sense  that  small-scale 
velocity  perturbations  are  not  considered.  Both  tomographic  and  scattering  studies  aim  at 
mapping  such  heterogeneities,  which  naturally  are  frequency  dependent.  Typical  "scatter¬ 
ing”  medium  paiameters  like  B-  and  S-velocilies  at  i  11/  for  the  crust  are  RMS  fluctua¬ 
tions  in  the  range  2-6  per  cent  and  correlation  distances  m  the  range  2-20  km  (Wu  and 
h'latte,  1988;  Charette,  1991).  Typical  correlation  functions  describing  such  media  are  of 
the  Gaussian,  exponential  or  the  von  Karman  type.  TTie  latter  can  describe  self-similar 
media,  by  which  is  meant  that  the  extent  of  velocity  perturbations,  calculated  over  equal 
logarithmic  intervals  of  wave  numbei,  remains  constant  over  a  range  of  scale  lengths 
(Frankel  and  Clayton,  1 986).  In  order  to  provide  a  better  understanding  of  Rg  propagation 
in  complex  media,  we  have  calculated  synthetic  seismograms  for  a  set  of  homogeneous 
and  inhomogeneous  models.  The  technique  used  is  that  described  by  Hestholm  et  al 
(1992'1  that  IS,  2-dimensional  (2D)  finite  difference  solutions  of  the  elastic  wave  equation. 
In  cases  of  inhomogeneous  upper  crust  media,  we  have  used  exclusively  2D  von  Karman 
functions.  As  mentioned,  it  is  considered  unrealistic  that  the  upper  crustal  low  velocity 
layer  (LVL)  should  have  a  uniform  thickness,  so  for  some  models  a  corrugated  interface 
was  introduced  in  the  form  of  a  1 D  von  Karman  function.  Details  on  the  various  models 
considered  are  given  in  Table  3.  Note  that  only  P-type  sources  are  used,  as  almost  all  Rg 
observations  stem  from  explosions  close  to  the  surface.  The  motivation  for  using  Rg  syn¬ 
thetics  was  as  follows: 

•  How  well  do  we  recover  model  parameters  when  inverting  the  synthetic  records  in 
the  same  manner  as  the  real  ones 

•  The  “response”  of  the  Rg  waves  to  model  heterogeneities. 

•  Rg  excitation  as  a  function  of  focal  depth;  this  problem  is  of  interest  in  the  context 
seismic  source  idcntilication  as  explosions  arc  confined  ai  most  to  1-2  km  depths. 
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Below  we  present  results  pertaining  to  the  above  class  of  problems  with  basis  in  syn¬ 
thetic  seismograms,  a  few  of  which  are  shown  in  Fig.  5  (Models  2,  3, 4)  and  Fig.  6  (Model 
2). 

Inversion  of  synthetic  seismograms 

From  synthetic  seismograms  as  exemplified  in  Fig.  5,  we  estimated  Rg  phase  veloci¬ 
ties  and  then  inverted  in  order  to  reuieve  the  original  model  parameters.  'I’he  outcome  of 
this  experiment  for  Model  2  (Table  3)  is  included  in  Table  2  for  the  distance  range  30-60 
km  and  depth  2  km.  The  differences  in  measured  synthetic  phase  velocities  compared  to 
theoretical  model  values  amounts  to  maximum  0.1  kms'^  At  shallow  depths  (1-2  km),  the 
Rg  fundamental  mode  would  be  dominant  over  the  higher  modes  for  entire  peiord  range 
considered  here  (e.g.,  see  Panza  et  al,  1973),  so  any  bias  effect  in  phase  velocity  measure¬ 
ments  is  likely  to  be  small.  Concerning  observational  data,  standard  practice  is  to  average 
phase  velocity  measurements  between  events  as  actually  done  here  and  also  in  Paper  1 . 

The  above  synthetic  experiment  has  been  instructive  as  it  demonstrated  that  Rg  phase 
velocity  inversion  provides  robust  estimates  of  top  crust  structural  parameters.  Also,  there 
is  obviously  a  trade-off  in  the  estimate  of  LVL  thickness  and  corresponding  shear  veloci¬ 
ties  at  least  for  observational  bandwidths  available  to  us. 

Rg  wavetrain  complexities 

Instructive  examples  are  shown  in  Fig.  5,  where  synthetics  in  two  different  passbands 
are  shown  for  Models  2,  3  and  4  in  Table  3.  At  relatively  low  frequencies  (0-2  Hz)  the  Rg 
wavetrain  is  almost  undeformed  even  when  companng  Model  2  and  4  synthetics.  How¬ 
ever,  at  higher  frequencies  (unfiltered  traces)  the  P-to-S  and  S-to-S  first  and  higher  order 
scattering  contributions  tend  to  mask  tlie  Rg  wavetrain.  Only  the  Airy-phase  arriving  atca 
25  sec  remains  clearly  discernible.  For  Model  5  (synthetics  not  shown)  these  effects  are 
even  more  pronounced.  The  Airy  phase  (period  ca  1  sec)  dominates  our  synthetics  which 
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is  typical  for  such  models.  Naturally,  epicenter  distance  ranges  are  of  importance,  as  for 
example  ASAR  at  450  km  (Table  1 )  only  records  the  long  period  part  of  the  Rg  (1 .7-4.5 
sec).  The  other  extreme  is  YKA  at  12-65  km,  where  Rg  periods  down  to  0.3  sec  are 
recorded.  In  view  of  our  synthetics,  this  would  imply  that  the  amount  of  scatters  in  top 
crust  near  YKA  must  be  moderate,  that  i.s,  must  be  less  than  the  RMS  of  4  per  cent  pre¬ 
sumed  for  Model  4. 

For  most  of  the  arrays,  the  source  distance  range  is  60-200  km  and  the  corresponding 
Rg  period  range  is  0.7- 1.7  sec.  For  distances  up  to  ca  100  km  interference  from  S-scatter- 
ing  wavelets  may  interfere  strongly  with  the  Rg  wavetrain.  The  main  effect  of  upper  crust 
heterogeneities  appears  to  be  to  attenuate  the  Rg  over  large  distances.  For  example,  the 
NORESS  array,  located  in  hilly  country,  hardly  records  Rg  beyond  100  km.  while 
ARCESS  records  Rg  out  to  ca  300  km.  From  Sweden  BSth  (1975,  1982)  reports  Rg  out  to 
600  km.  Very  efficient  topography  related  P-to-Rg  scattering  has  been  reported  by  Bannis¬ 
ter  et  al  (1990).  Our  synthetics  were  computed  without  incorporating  intrinsic  attenuation 
in  the  models  and  are  not  aimed  at  addressing  this  problem.  On  the  other  hand  they  illus¬ 
trate  that  complex  media  (Models  3  and  4)  may  give  rise  to  mode  conversion  and  other 
scattering  effects. 

Rg  excitation  as  a  function  of  depth. 

The  results  in  Fig.  5  were  obtamed  for  a  fixed  source  depth  at  2  km.  Moving  the 
source  depths  to  4  and  6  km,  respectively,  had  rather  profound  effects  on  the  Rg  excitation 
efficiency  as  shown  in  Fig.  6.  Essentially,  for  depths  below  2  km  very,  very  little  Rg  exci¬ 
tation  takes  place  in  the  0.5  -  2.0  see  signal  band.  In  other  words,  Rg  observations  imply 
that  the  source  must  be  very  shallow.  The  array  data  used  m  this  study  stem  from  chemical 
explosions  near  the  surface,  except  those  from  the  Alice  Springs  array.  In  the  latter  case, 
earthquake  rupturing  of  the  free  surface  took  place. 
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4  Discussion 


In  this  study  we  have  used  Rg  recordings  from  seven  arrays  located  on  four  continents 
to  explore  structural  inhomogeneities  in  the  uppermost  crust  in  their  respective  siting 
areas.  The  common  result  feature  is  the  existence  of  a  low-velocity  layer  (LVL)  in  the 
uppermost  part  of  the  crust.  This  may  to  some  extent  be  an  artefact  of  the  inversion  models 
used;  we  have  not  considered  models  with  gradual  velocity  increases  in  layer(s)  over  a 
halfspace.  Anyway,  the  shear-wave  estimates  for  the  halfspace  are  in  average  3.55  kms'* 
and  remarkably  consistent  between  the  widely  separated  arrays  (details  in  Table  2).  EKA 
is  exceptional  in  this  regard  (3.10  kms'*),  which  is  attributed  to  its  siting  on  sedimentary 
rocks.  The  shear  velocity  within  the  presumed  LVL  is  reasonably  consistent  (Table  2)  and 
in  average  2.87  kms  '.  Estimated  layer  thicknesses  vary  considerably,  being  thin  for  YKA 
(0.12  km)  and  ARAO  (0.45)  and  very  thick  for  ASAR  (1.6  km).  We  have  also  estimated 
shear  velocities  using  a  model  with  fixed  layer  thicknesses  over  a  halfspace.  The  results 
here  (Table  2)  demonstrate  that  LVL  thickness(es)  and  corresponding  shear  velocities  are 
not  well  resolved. 

The  observational  data,  that  is,  Rg  phase  velocity  dispersion  curves,  vary  somewhat 
from  one  event  to  another,  hence  averaging  over  a  few  events  was  performed.  The  reason 
for  this  is  attributed  to  interferences  from  S-scattering  wavelets.  This  has  been  visualized 
through  synthetic  seismogram  analysis  using  a  suite  of  homogeneous  and  inhomogeneous 
structural  models  (Table  3).  It  has  not  been  possible  to  mimic  extreme  observational  fea¬ 
tures  like  the  vanishing  of  Rg  on  one  sensor  while  another  one  500  to  1000  m  away  exhib¬ 
its  clear  Rg  phases.  Also,  an  LVL  model  would  imply  the  existence  of  a  clear  Airy  phase 
with  typical  features  like  the  synthetics  displayed  in  Fig.  5.  Observationally,  the  Rg  wave¬ 
form  IS  more  complex  (seldom  very  clear  Airy  phases),  which  may  reflect  the  effects  a 
gradual  velocity  increase  in  the  uppermost  crust  or  interferences  from  S-scattering  wave¬ 
lets. 


78 


Hie  main  result  of  tins  study  is  that  a  low-velocity  layer  appears  to  exist  in  the  upper¬ 
most  part  of  the  continental  crust  and  that  this  is  a  global  phenomenon.  A  common  physi¬ 
cal  explanation  is  that  of  a  relative  abundance  of  cracks  near  the  surface  as  well  as  a  water 
content  that  is  relatively  high  in  the  uppermost  1-2  km  of  the  crust  (e.g.,  see  Fritz  and 
Frape,  1987).  This  provides  a  mechanism  for  lowering  .seismic  velocities.  Results  of  a  lab¬ 
oratory  study  of  such  problems  have  been  reported  by  Stesky  (1985),  who  found  that  frac¬ 
tures  in  rocks  under  confining  pressures  up  to  200  MPa  decreased  P-  and  S-velocities  on 
the  order  of  10-20%.  The  effect  was  greater  at  lower  pressures,  with  higher  averaged  num¬ 
bers  of  fractures,  and  in  rocks  containing  a  lower  microcrack  porosity.  Although  peuolog- 
ical  factors  are  claimed  to  be  important  on  the  basis  of  Rg  studies  in  mainly  sedimentary 
rock  environments  (MacBeth  and  Bunon,  1986,  1987;  Saikia  et  al,  1990,  among  others), 
this  is  not  obvious  from  our  results  stemming  from  Rg  observations  in  crystalline  rock 
environments  (except  EKA). 

Finally,  we  want  to  comment  on  the  practical  aspects  of  our  results.  Clearly,  a  pro¬ 
nounced  LVL  on  top  of  the  crust  would  affect  the  waveform  of  high-frequency  P-  and  S- 
waves  propagating  through  such  a  layer.  For  example,  Lokshtanov  et  al  (1991)  demon¬ 
strated  that  more  consistent  slowness  estimates  for  P  were  obtained  if  the  LVL  effect  was 
taken  into  account.  S-waves  having  relatively  short  wavelengths  are  likely  to  be  affected 
even  more  by  such  a  low-velocity  layer,  which  observational ly  is  manifested  by  rather 
complex  particle  motion  patterns  (Roberts  and  Christoffersson,  1990).The  observational 
existence  of  Rg  waves  implies  that  the  source  depth  should  be  less  than  say  2-3  km  (see 
Fig.  6).  Thus  in  many  contexts,  Rg  observations  can  be  used  to  discriminate  between  man¬ 
made  explosions  and  real  earthquakes  (Kafka,  1990).  In  practice,  this  would  imply  that 
more  effective  means  for  extracting  Rg  parameters  from  the  seismic  records  need  to  be 
developed. 
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5  Concluding  remarks 


Using  Rg  recordings  from  seven  arrays  on  four  continents,  we  have  demonstrated  that 
an  uppermost  crust  low-velocity  layer  appears  to  be  a  global  continental  phenomenon.  The 
halfspace  average  shear  velocity  of  3.55  kins''  is  remarkably  consistent  among  the  arrays 
and  the  same  applies  roughly  to  the  LVl,  shear  velocity  average  oi  2.87  kms  '.  However, 
the  LVL  thicknesses  vary  considerably,  which  to  some  extent  may  reflect  poor  resolution 
for  this  parameter  in  the  observational  data.  Through  syntheuc  seismogram  analysis,  we 
have  demonstrated  that  structural  inhomogeneities,  modelled  as  a  (self-similar)  von  Kar- 
man  medium,  could  to  some  extent  mask  the  Rg  wavetrain  through  shear  wave  scattering 
contributions.  For  events  with  focal  depths  exceeding  2-3  km,  Rg  excitation  would  be 
minimal,  and  thus  Rg  ooservations  have  the  potential  of  being  a  powerful  earthquake/ 
explosion  discriminant  for  epicentral  di.stances  of  a  few  hundred  kilometers.  The  existence 
of  the  l.VL  naturally  contiibuios  to  the  complexity  of  local  seismogram  records  where 
high  signal  frequencies  (>4  Hz)  are  dominant. 
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Array 

Lat 

(deg) 

Long. 

(deg) 

Geol.  Siting 

Source 

Ranges 
(km)  (deg) 

Rg-Periods 

(sec) 

NORESS-NRAO 
(S.  Norway) 

60.74N 

11.54E 

Piecanibrian 

Granite-gneiss 

1  EX 

34 

314 

0.7- 1.5 

ARCESS-ARAO 
(N.  Norway) 

69.53N 

25.5  IE 

Precambnan 

Gabbro 

2  EX 

170-190 

83-96 

0.65-1.70 

GERESS-GEC2 
(Bavaria,  Germany) 

48.84N 

13.70E 

Precambnan 

Granite-gneiss 

4  EX 

70-90 

130-325 

0.6-1. 4 

Eskdalemuir-EKA 

(Scotland) 

55.38N 

3  13W 

Silunan  shales 

3  EX 

60 

279 

0.6-21 

Garibidinaur-GBA 

(India) 

13.60N 

77.44E 

Archean  gneiss 

3  EX 

110-215 

125-330 

0.8-2.0 

Yellowknife-YKA 
(NW  Canada) 

62.49N 

114.61W 

Archean  gneiss 

2  EX 

12-65 

90-300 

0.3-2.0 

Alice  Spnngs-ASAR 

23  67S 

133  90E 

Piecanibrian 

2EQ 

450 

6-7 

1. 7-4.5 

(C.  Australia) 


Table  1.  Location  of  the  arrays  used  in  this  study.  Brief  information  on  geological  siting 
and  the  events  analyzed  are  given  (EX  =  explosion;  EQ  =  earthquake).  Distance 
ranges  of  the  events  are  given  in  km  and  azimuth  range  in  degrees. 
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Ono  layer  model 

'Fwo  layers  model 

Array 

/?! 

/ii 

02 

a 

ftx 

02 

06 

a 

(km/s) 

(km) 

(km/s) 

(kin/s) 

(kin/s) 

(km/s) 

NRAO 

3.02 

1.07 

3.52 

0.0014 

2.80 

3.36 

3.59 

0.0003 

ARAO 

2.69 

0.45 

3.65 

0.0006 

3.02 

3.50 

3.62 

0.0019 

GEC2 

2.93 

1.02 

3.54 

0.0012 

2.59 

3.35 

3.63 

0.0018 

EKA 

2.59 

0.71 

3.10 

0.0028 

2.41 

3.09 

3.12 

0.0053 

GBA 

2.86 

0.76 

3.51 

0.0238 

2.57 

3.52 

3.52 

0.0364 

YKA 

2.81 

0.12 

3.44 

0.0005 

3.30 

3.38 

3.42 

0.0011 

ASAR 

2.96 

1.60 

3.68 

0.0275 

2.91 

2.91 

3.70 

0.0320 

SYN 

2.68 

1.23 

3.43 

0.0062 

2.51 

2.88 

3.48 

0.0382 

Table  2.  Results  from  inversion  of  dispersion  data.  P  is  S-wave  velocities  and  h  is  layer 
thickness.  For  Model  2  the  layer  thicknesses  were  fixed  to  hj  =  0.5  km  and  h2  = 

n 

1 .0  km.  a  is  a  data  misfit  function  defined  as  o  =  1000  ( v“‘  -  m"‘)  ^  where  v  and 

I  =  1 

u  are  observed  and  calculated  phase  velocities,  respectively  and  n  is  the  number  of 
observations.  In  the  last  line  (array  =  SYN)  results  of  analysis  and  inversion  of  the 
synthetic  data  are  shown. 


Synth. 

model 

Model 

type 

Model 

dist. 

So 

type 

iirce 

depth 

(km) 

Ba 

01 

(km/s) 

isic  mot 
hi 

(km) 

iel 

02 

(km/s) 

von 

u 

Karmai 

a 

(km) 

n  medium 

rms 

1 

h. space 

no 

P 

2/4/6 

- 

- 

- 

- 

- 

2 

LVL 

no 

P 

2/4/6 

2.82 

1.4 

3.55 

- 

- 

- 

3 

LVL 

Cl 

mm 

2/4/6 

2.82 

1.4 

3.55 

0.5 

10 

0.1km 

4 

LVL 

RM 

p 

2/4/6 

2.82 

1.4 

3.55 

0.0 

5 

4% 

5 

LVL 

Ani. 

p 

2/4/6 

2.82 

1.4 

3.55 

0.0 

5/1 

4% 

Table  3.  Summary  of  the  models  used  for  computation  of  finite  difference  synthetics.  In 
the  column  “Model  disturbance”  Cl  means  corrugated  interface  generated  by  1 D 
randomization  of  layer  thickness  using  a  von  Karman  correlation  function  with  indi¬ 
cated  parameters  (order,  correlation  distance,  and  rms  variation).  RM  means  random 
media  (velocity  vanations)  in  the  LVL  with  given  correlation  function.  “Am.”  is  for 
a  model  with  “apparent  anisotropi”,  that  is,different  correlation  distances  in  the  hor¬ 
izontal  and  vertical  direction.  (Results  of  this  last  model  is  not  shown  here). 
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Fig.  1.  Typical  waveforms  fora  shallow,  local  (distance  about  215  km)  event  showing 
clear  Rg  waves.  The  seismograms  are  from  seven  stations  along  the  “red”  line  of  the 
Garibidinaur  anay  in  India.  The  Rg  phases  starting  at  about  70  seconds  have  the 
highest  amplitudes  in  the  seismograms  and  aie  chiuacterized  by  low  frequencies  and 
a  dispersive  wavetrain  (with  gradually  increasing  frequency).  Also  seen  is  the  S- 
phase  at  about  55  sec  and  the  P-phase  at  30  sec. 
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Fig.  2.  Dispersion  analysis  of  the  event  shown  in  Fig.  1  (all  20  channels  are  used).  The 
solid  line  is  the  dispersion  curve  used  for  inversion. 
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PERIOD  (S) 

Fig.  3.  The  averaged  dispersion  curves  extracted  for  the  seven  arrays  used  in  this  study. 
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Fig.  4.  S-velocities  versus  depth  for  the  two  types  of  models  considered.  On  the  left  side 
results  from  Model  1  inversion  are  shown  fl  layer,  variable  velocities  and  thick¬ 
ness),  while  Model  2  results  (2  layers,  variable  velocities,  fixed  thicknesses)  are 
shown  to  the  right. 
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Fig.  5.  Synthetic  seismograms  computed  by  2D  finite  difference  elastic  wave  modelling 
for  different  models.  In  all  three  cases  the  explosion  source  is  located  at  60  km  dis¬ 
tance  and  2  km  depth.  The  trace  labelled  LVL_1  was  produced  with  a  flat  homoge¬ 
neous  low  velocity  layer.  LVL_2  is  for  a  model  with  random  variations  in  layer 
thickness,  while  LVL_3  is  for  a  fiat  layer  with  random  variations  in  layer  velocity. 
Further  model  details  are  given  in  Table  3.  Below  each  of  the  three  traces,  a  low- 
pass  filtered  (0-2  Hz)  version  is  shown.  Note  that  the  low  frequency  Rg  waves  are 
almost  unaffected  by  the  random  variations. 
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Fig.  6.  Synthetic  seismograms  for  different  source  depths  (explosion  source,  homoge¬ 
neous  LVL  model  and  60  km  distance  for  all  cases).  Source  depths  are  2, 4  and  6 
km.  Rg  excitation  is  seen  to  diminish  very  rapidly  with  depth. 
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Abstract 

In  this  study  we  report  on  efforts  in  using  synthetic  seismograms  for  a  better  under¬ 
standing  of  wave  propagation  in  complex  media.  Firstly,  a  comprehensive  description  is 
given  of  the  numerical  finite  difference  (FD)  technique  used  to  solve  the  elastic  wave 
equation.  Our  FD  schemes  have  been  tested  on  2-dimensional  (2D)  crust/lithosphere  mod¬ 
els  of  varying  complexities.  The  first  model  considered  was  a  simple  one-layered  homoge¬ 
neous  crust,  then  the  model  was  made  more  complex  by  adding  corrugations  to  the  Moho- 
interface.  Further  model  inhomogeneities  introduced  were  randomized  velocities  with  cor¬ 
relation  functions  of  the  von  Karman  type.  The  main  result,  when  comparing  to  observa¬ 
tional  records  at  local  distances,  is  that  c.  4  per  cent  RMS  velocity  fluctuations  in  the  crust 
IS  required  to  ensure  sufficient  strong  and  persistent  coda  exiiation.  However,  coda  coher¬ 
ency  remains  relatively  high  which  we  attribute  to  the  limitation  of  using  2D-schemes  and 
also  that  scattering  by  free  surface  topography  is  not  taken  into  account.  The  relatively 
small  amplitudes  of  the  Pn-phase  con.pared  to  later  P-phases  in  our  synthetics  may  reflect 
a  lack  of  sub-Moho  velocity  gradients  in  the  model.  Finally,  S-to-P  conversions  appear  to 
be  more  efficient  than  P-ti'-S,  which  may  prove  to  be  a  diagnostic  feature  for  better  event 
classification. 
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1  Introduction 


The  understanding  of  high-frequency  (  >  1  Hz)  seismic  wave  propagation  in  the  crust/ 
lithosphere  structural  system  remains  problematic  to  seismologists  and  exploration  geo¬ 
physicists  alike.  The  observational  manifestations  of  this  enigma  are  easy  to  illustrate;  at 
local  distance  ranges  both  P  and  S  waves  are  accompanied  by  prolonged  coda  waves 
which  are  not  accounted  for  by  standard  structural  models,  likewise,  array  recordings  of 
teleseismic  events  exhibit  intrinsic  P-wave  time  and  amplitude  anomalies,  also  between 
sensors  spaced  only  1-2  km  apart  (Haddon  and  Husebye,  1978).  Howevei,  complex  lay¬ 
ered  models  in  combination  with  wave  propagation  schemes  tied  to  ray  tracing,  Gaussian 
beams,  WKBJ,  perturbation  methods,  etc.,  cannot  realistically  model  observational  fea¬ 
tures  of  the  above  kind.  The  rapid  P-wave  time  and  amplitude  fluctuations  across  the  large 
arrays  NORSAR  (SE  Norway)  and  LASA  (Montana)  are  reminiscent  of  scattering.  The 
first  successful  attempts  at  seismological  modelling  of  such  phenomena,  using  the  first 
order  (acoustic)  scattering  theory  of  Chernov  (1960),  were  Aki  (1973)  (see  also  Bertheus- 
sen  et  al,  1975;  Platte  and  Wu,  1988,  among  others).  In  this  context,  the  scattering  medium 
IS  described  statistically  in  terms  of  first  and  second  order  moments  and  the  correlation 
distance  length.  Naturally,  in  the  solid  earth  scatterers  must  have  deterministic  locations, 
as  elegantly  dcmonstiatcd  by  Aki  ct  al  (1977),  using  tlieii  novel  tomographic  technique  in 
analysis  of  array  P-time  residuals. 

Effons  to  exploit  the  information  potential  of  the  coda  waves  from  small,  local  earth¬ 
quakes  have  met  with  relatively  little  success  in  terms  of  locating  specific  scattering 
sources.  On  the  other  hand,  Aki  and  Chouet  (1975)  demonstrated  that  characteristic  coda 
features  like  duration,  decay  rates  and  scattering  attenuation  can  be  modelled  in  terms  of  a 
single  backscattering  theory  (weakly  inhomogeneous  media  -  the  Born  approximation 
valid).  Further  elaborations  on  the  single  scatter  theory  come  from  the  work  of  Wu  and 
Aki  (1985)  and  Wu  (1985.  1988).  who  split  the  scattering  contnhuiions  in  two  parts, 
namely,  a  back-scattering  tcim  stemming  from  impedance  perturbations  and  a  forward- 
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scattering  term  stemming  from  velocity  perturbations.  Dainty  and  Toksoz  (1977)  consid¬ 
ered  similar  problems  but  using  diffusion  theory. 

The  basic  problem  in  most  studies  of  wave  propagation  in  complex  media  is  that  scat¬ 
tering  theory  incorporating  the  Born  approximation  is  not  necessarily  valid,  while  conve¬ 
nient  analytical  solutions  of  the  elastodynamic  wave  equation  are  seldom  at  hand.  In  this 
context  it  is  naturally  to  aim  for  the  numerical  solutions  of  the  wave  equation  which  have 
been  around  for  some  time  (Alterman  and  Kaial,  1968).  Finite  difference  solutions  here 
can  handle  complex  media  and  at  the  same  time  properly  account  for  multiple  scattering 
effects,  although  severe  computer  requirements  still  limit  the  practical  usage  of  such  pow¬ 
erful  modelling  tools.  Most  of  the  works  published  so  far  appear  to  be  aimed  at  testing  the 
validity  of  the  numerical  codes  used,  exploring  the  practical  applicability  of  the  Born 
approximation  in  scattering  studies  and  recently  investigating  the  range  of  “scattering” 
functions  appropriate  for  crust/lithosphere  medium  descriptions.  Works  of  relevance  here 
are  Macaskill  and  Ewart,  1984;  Frankel  and  Clayton,  1986;  McLaughlin  and  Anderson, 
1987;  Dougherty  and  Stephens,  1988;  Toksoz  et  al,  1990. 

In  tins  study  we  concentrate  on  generating  FI)  synilietics  with  basis  in  the  clastody- 
namic  wave  equation,  and  limit  ourselves  to  2D  crust/lithosphere  models.  The  starting 
point  is  details  on  the  velocity-stress  formulation  used  in  the  numerical  solution  of  the 
wave  equation.  Then  source  functions  incorporated  and  scattering  medium  functions  con¬ 
sidered  would  be  described.  In  the  latter  case  it  is  of  some  importance  whether  medium 
properties  are  non-uniform  within  the  lithosphere  as  claimed  by  Flatte  and  Wu  (1988). 
Results  are  in  terms  of  synthetics  for  realistic  crustal  models,  which  are  subsequently  dis¬ 
cussed  in  terms  of  characteristic  features  of  NORESS  array  recordings. 
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2  Discretization  of  the  wave  equation 


The  basic  equations  governing  wave  propagation  in  a  continuous  elastic  medium  are 
the  momentum  conservation  and  the  stress-strain  relation.  Following  Achenbach  (1984), 
in  the  velocity-stress  formulation,  these  are  given  by 


P3,'' 


s}  ~  'j  r.'-i . J 

|j  Oy,  =  v,h-52^  vp  y./=l . J  J*l 


(1) 


(2) 


(3) 


where  Einstein’s  summation  convention  is  used.  7  is  the  dimensionality  of  the  problem,  p 
is  density,  and  X  and  p  are  Lame’s  parameters.^  are  body  forces  and  v^  and  Ojj  are  parti¬ 
cle  velocities  and  stresses,  respectively. 


Spatial  partial  differentiation  is  achieved  through  cost-optimized,  dispersion-bounded, 
high-order  finite  difference  operators  on  a  staggered  grid.  For  time  stepping  a  leap-frog 
technique  is  used.  The  discretization  of  the  elastodynamic  equations  with  two  staggered 
numerical  space  differentiators,  6*  ,  applied  as  in  Levander  (1988)  to  stresses  and  particle 
velocities  leads  to: 


p)  [V*  (t  +  M/2)  -V^  (t-M/2)]  =  At{Fj  (l)  -i-bjSjjil)  +  (0} 

j,l=l . J 

J 

S,:(r  +  Ar)  -S.  (0  =  XAt  V  b'  (l  +  At/2)  +  2\iAib  V*  +  M/2) 

JJ  JJ  ^  '  J  J 

r  =  1 

jj=\ . J 

s;;(/+Af)-5;;(o  =\ij;At{bjv^{t)+bivjit)}  j,i=:  i . j-j*i 
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with 


/;  (/)  =  V . (X  +  h/2.  /)  F]  (/)  =  /^.  (X  +  h/2.  t) 


Sjj ( t)  =  Ojjix,  l)  Sj^  ( t)  =  Ojjix  +  h y 2  +  h/2,  0 


Pj 


=  p(x  +  h/2)  X,  =  X(x)  [I  =  |l(x) 


and 


ji 


=  p(x  +  h./2  +  h/2) 


^-'^2  ^(x  +  /h.)  -i?(x- (/- l)h  )  dq 

5;,(x)  =  - - 

/  =  1  ^  > 


<7(x+ (/- l)h.) -<7(x-/h  )  dq 

5-«(x)  =  - ± - |(x-h/2) 

t  =  1  J  J 

Here  hy  is  the  unit  vector  m  the  jth  direction,  X,  n  and  are  defined  at  the  nodes  of  the 
Cartesian  mesh,  pj,  V*  and  F*  are  defined  at  the  links  connecting  the  nodes  and  5/  and 
g/  are  defined  at  the  centers  of  the  “obliquities’.;  6-  are  numerical  differentiators  of 
coefficients  .q  is  velocity  or  stress  and  L*  is  the  length  of  the  operator. 


For  the  numerical  dispersion  relations,  the  stability  limit  and  bandwidth  introduced  by 
the  discretization,  the  reader  is  referred  to  Sguazzero  et  al  (1989).  The  computational  cost 
(in  floating  point  operations)  of  an  elastic  numerical  simulation  is  given  by 


c„=  y(%  +  <  +  2)+7J-ll(N„/iJ^+'a;>  (4) 

Here  ^  2  is  the  number  of  grid  points  per  shortest  wavelength  needed  to  model  the 
wavefield,  =  3  (L*  /2)  -  1  is  the  number  of  floating  point  operations  needed  to  com¬ 
pute  the  numerical  derivative  and  is  the  number  of  shortest  wavelengths  in  each  direc¬ 
tion  of  the  computational  domain.  7  =  1 , 2  or  3  is  the  dimensionality  of  the  problem, 

=  cht/dix  is  the  Courant  number  of  this  explicit  scheme  and  c  is  the  relevant  propaga¬ 
tion  speed. 
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spatial  partial  differentiation  is  performed  after  the  optimization  of  (4)  with  respect  to 
,  subject  to  various  constraints  (for  details,  see  Sguazzero  et  al,  1989).  This  process 
leads  to  the  optimized  operators  utilized  in  our  code. 

Absorbing  Boundary  Conditions 

By  necessity,  the  numerical  modeling  limits  the  medium.  To  reduce  artificial  reflec¬ 
tions  from  the  numerical  boundaries,  the  velocities  and  stresses  are  multiplied  by  expo¬ 
nentially  decreasing  terms  near  the  edges.  Using  this  procedure,  the  reflections  are  still 
visible  to  a  certain  extent  by  the  time  the  wave  has  been  reflected  and  propagated  back  into 
the  computational  domain.  This  causes  errors.  Therefore  this  treatment  of  boundary  reflec¬ 
tions  requires  a  relatively  long  distance  in  each  spatial  direction,  and  thus  in  3D  becomes 
computationally  intensive.  A  typical  application  demands  a  run  time  storage  of  700-800 
MBytes. 

A  recently  developed  alternative  to  this  procedure  has  been  forwarded  by  Higdon 
(1990,  1991).  He  introduced  an  operator  to  minimize  the  reflections  at  the  computational 
boundary.  At  a-  =  0  this  reads 


n 


(5) 


;=  1 

which  will  absorb  perfectly  a  plane  wave  traveling  towards  the  boundary  at  angle  a.  and 
speed  Cj.  m  reflects  the  accuracy  of  the  operator,  and  is  the  number  of  angles  for  which 
absorption  is  perfect.  Similar  operators  are  used  on  the  other  boundaries.  A  version  of  this 
algorithm  has  been  successfully  applied  to  various  body  geometries  and  structures.  With  a 
suitable  positioning  of  source  and  receivers,  this  procedure  could  halve  the  distance  neces¬ 
sary  in  each  spatial  direction.  The  required  3D  run  time  storage  for  the  computational 
domain  is  then  reduced  to  1/8  of  the  above.  Nevertheless,  it  is  the  method  of  exponentially 
decreasing  terms  for  damping  “edge  effects”  that  has  been  used  in  our  2D  synthetics.  The 
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accuracy  of  this  method  proved  to  be  better  than  that  of  Higdon  (1990).  In  the  latter  case, 
the  problem  appears  to  be  that  the  number  of  incident  angles  for  which  the  absorption  is 
good  is  limited. 

Free  Surface  Boundary  Conditions:  On  the  top  surface,  we  use  the  vanishing 
.stress  conditions  for  a  free  boundary 


n  •  T  -  0  (6) 

Here  h  is  the  outward  normal  unit  vector  of  the  surface  and  T  is  the  stress  tensor.  To  get 
computationally  tractable  conditions,  we  assume  the  free  top  surface  to  be  locally  plane, 
with  «  =  i-,  where  k  is  the  unit  vector  in  the  vertical  z-direction.  x  and  y  are  horizontal 
coordinates.  Eq.  (6)  then  leads  to 


a 


zx 


(7) 


Across  the  discontinuity  at  the  top  surface  the  momentum  conservation  equation  (1)  is 
not  valid.  Eq.  (7)  put  into  the  versions  of  eqs.  (2)  and  (3)  for  the  stresses  o  and 
leads  to  the  system 


dv,  dv^ 

dz  Bx 


(8) 


dz  dy 


(9) 


d^, 

dz 


A,  +  2p.  dx  dy 


(10) 


which  is  closed  with  respect  to  the  three  velocity  components.  Once  these  are  solved  for,  it 
is  possible  to  find  the  remaining  stresses  and  on  the  surface  from  the  appropri¬ 

ate  versions  of  eqs.  (2)  and  (3). 
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Note,  the  above  boundary  conditions  arc  strictly  valid  for  a  locally  plane  free  surface  and 
not  for  one  exhibiting  topographic  relief.  The  latter  case  is  rather  problematic  even  for  2D 
models  (Jih  et  al,  1988),  and  at  present  such  an  option  has  not  been  incorporated  in  our 
numerical  scheme. 

Discretization  of  Free  Surface  Conditions:  The  dependent  variables,  Lame’s 
parameters,  the  density  and  volume  forces  are  defined  on  a  staggered  grid  in  accordance  to 
Levander  (1988).  Because  of  the  way  in  which  the  particle  velocities  and  stresses  are 
defined  on  this  grid,  it  turns  out  that  eqs.  (8)-(10)  can  be  resolved  as  an  explicit  system  of 
equations  once  the  interior  points  have  been  solved  for.  If  only  velocities  are  needed,  this 
suffices  as  the  system  is  closed.  Otherwise  eqs.  (2)  and  (3)  can  subsequently  be  used  to 
calculate  the  remaining  free-surface  stresses. 

In  our  application  a  cost-optimized,  dispersion-bounded  8th  order  operator  of  length  8 
IS  chosen  as  spatial  differentiator.  It  corresponds  to  a  sampling  frequency  of  3  per  mini¬ 
mum  wavelength,  and  it  has  a  relative  error  bound  of  1 .9%  for  the  group  velocity.  Because 
of  the  length  of  the  operator,  other  devices  have  to  be  used  on  the  layers  near  all  bound¬ 
aries.  Second  order  staggered  finite  differences  given  earlier  as  bj  with  L-  =  2  are  there¬ 
fore  employed  at  all  4  outermost  points  of  the  domain.  These  spatial  differentiators  are 
also  used  for  discretizing  the  system  (8)-(10). 

Recently,  our  scheme  has  been  improved  by  increasing  the  accuracy  near  the  free  sur¬ 
face.  Here  we  use  second-order,  single-sided,  staggered  finite  difference  operators  for  the 
discretization  of  eqs.  (8)-(10).  Below  the  surface  the  operator  order  gradually  increases 
with  depth,  that  is,  from  second  to  fourth  and  sixth  order,  and  these  central,  uniform  oper¬ 
ators  are  employed  on  a  staggered  grid  (Fornberg,  1988).  In  this  way,  a  high  level  of  accu¬ 
racy  is  maintained  at  the  surface.  For  the  Courant  number  used,  these  operators  are  stable 
even  when  simulating  very  short  wavelength  Rg  phases.  These  are  fundamental  mode 
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Rayleigh  waves  propagating  along  the  surface  with  phase  velocities  in  the  range  2.5- 
3.5  kms'^  for  periods  0. 5-4.0  sec  (e.g.,  see  Ruud  et  al,  1992). 

3  S^furce  types  used  in  seismic  modeling 

Below  we  give  the  body  force  equivalents  for  the  two  source  types  used  in  this  study. 
We  assume  that  the  source  extent  is  small  compared  to  the  dominating  wave  lengths  so 
that  we  can  use  a  point  source  approximation.  Here  s{t)  is  the  source  time  function  and 
6(x)  is  the  Dirac  delta  function. 

Explosion  Source  (P- source) 

/,(.r,z,/)  =  ^^b(x-XQ)b(2-ZQ)s(t)  (11) 

f.(x,z,t)  =  -b(x~X(i)  ^J)iz-ZQ)s(t)  (12) 

This  force  field  has  zero  curl  ( V  x  f  =  0)  so  that  no  S-waves  are  generated  at  the  source. 

Rotation  Source  (S-source): 


/,  (a,  2, /)  =  +b{x-XQ)  ^-b{z-ZQ)  s  (t)  (13) 

f^{x,z,t)  =  -^-bix-XQ)biz-ZQ)s{t)  (14) 

This  force  field  has  zero  divergence  (V  f  =  0)  so  that  no  P-waves  are  generated  at  the 
source.  Unlike  the  explosion  source,  this  is  not  an  internal  source  because  the  angular 
moment  of  the  force  field  is  different  from  zero. 

In  practice  an  approximation  must  be  used  for  the  delta  function.  We  have  used  a 
Gaussian  function  (with  a  =:  2\h^■  =  2\h)\ ); 
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(15) 


6(x)  =  ,  '  e 
j2no 


k 

2a^ 


,6  (A) 


a 


For  s(t)  we  have  used  a  function  proportional  to  ^  6  (/)  with  o  •«  0.036  s. 

of 


4  Scattering  representation  of  the  lithosphere 

Seismic  wave  scattering  is  a  complex  phenomenon  which  depends  on  the  size,  distri¬ 
bution  and  magnitude  of  the  heterogeneities  in  the  lithosphere  which  in  general  are 
unknown.  Traditionally,  the  earth  is  often  modelled  as  a  simple  stratified  medium,  each  of 
the  strata  having  constant  physical  properties.  Seismograms  from  these  models  tend  to 
match  the  gross  features  of  observational  records  but  lack  the  variations  in  amplitude  and 
travel  time  and  coda  waves  accompanying  major  phase  arrivals.  These  features  are  symp¬ 
tomatic  of  scattering  from  small-scale  changes  in  velocity  and/or  density. 

In  the  scattering  literature,  heterogeneous  media  are  commonly  described  in  terms  of  a 
few  physical  parameters  like  thickness  of  the  scattering  layer,  heterogeneity  correlation 
distance  a  (in  case  of  anisotropy  a^,  a^)  and  heterogeneity  fluctuation  (RMS  variation).  We 
may  limit  fluctuations  to  either  velocity  or  density  and  also  introduce  corrugated  layer 
boundaries.  Heterogeneous  media  realizations  are  often  represented  by  random  fields 
where  complexities  can  be  expressed  in  terms  of  few  low  order  statistical  moments  based 
on  the  above  scatter  parameters  a  and  RMS.  In  this  context,  the  Gaussian,  the  exponential 
and  the  von  Karman  correlation  functions  have  become  popular  among  seismologists.  The 
fall-off  rate  of  the  spectra  controls  the  amount  of  roughness  in  the  realizations.  The  von 
Karman  function  appears  most  suitable  for  describing  lithospheric  heterogeneities  (e.g., 
Frankel  and  Clayton,  1 986;  Toksoz  et  al,  1988;  Charette,  1991).  Note,  we  may  have  sev¬ 
eral  medium  realizations  with  basis  in  the  von  Karman  function  reflecting  a  particular 
choice  of  the  v  parameter  in  Table  1 .  For  example  for  v  =  0.5  it  simplifies  to  an  exponen¬ 
tial  function  while  for  v  =  0.0  the  medium  is  characterized  with  heterogeneities  that  are 
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self-similar  for  ka  >  1  (k  is  wave  number).  With  “self-similar”  is  meant  that  the  standard 
deviation  of  the  medium,  calculated  over  equal  logarithmic  intervals  of  wave  number 
remains  constant  over  a  range  of  scale  lengths  (hrankel  and  Clayton,  1986).  Note  that 
Plane  and  Wu  (1988)  in  their  3D  scattering  analysis  of  NORSAR  P-waves  used  a  band- 
limited  power  law  function,  which  is  quite  similar  to  that  of  the  exponential  function. 

In  the  present  work  we  have  exclusively  experimented  with  von  Karman  media  with 
V  =  0.0  and  v  =  0.5.  Likewise,  we  have  used  a  similar  1 D  von  Karman  function  to 
describe  a  corrugated  Moho.  In  general,  the  scattering  medium  is  isotropic,  but  can  easily 
be  made  “apparently”  anisotropic  by  simply  using  different  correlation  lengths  along  the 
X-  and  z-  axes,  respectively. 


5  Synthesizing  wave  propagation  in  the  crust  and  lithosphere 

The  problem  at  hand  is  perhaps  most  easily  illustrated  through  display  of  real  seismic 
records  as  in  Figs.  1  and  2.  The  first  one  shows  3-component  NORESS  recordings  from  a 
ripple-fired  quarry  blast  300  km  away  to  the  southwest.  The  surprising  feature  here  is  that 
the  dominant  wavefield  amplitude  is  associated  with  shear  waves  (Lg  phases)  on  the  trans¬ 
verse  compionent!  In  Fig.  2  a  refraction  profiling  section  from  EUGENO-S  (Gregersen, 
1991)  is  shown.  The  outstanding  feature  here  is  that  the  coda  dominates  the  wavefield 
after  the  PmP-phase  arrival  (PmP  -  P-wave  reflected  from  the  Moho).  To  us  it  is  rather 
speculative  to  attempt  to  identify  specific  arrivals  within  the  P-coda.  In  the  following,  we 
present  a  succession  of  crust/lithosphere  synthetics  for  a  suite  of  models  of  increasing 
structural  complexities.  The  basic  model  parameters  are  a  crustal  thickness  of  30  km,  and 
P-velocities  of  respectively  6.15  kms  ’  and  8.15  kms’’  (sub-Moho).  P-  and  S-velocities 
are  related  through  a  Poisson’s  ratio  of  0,25. 
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Homogeneous  crust  with  different  velocity  distributions 

The  homogeneous  crustal  model  considered  is  shown  in  Fig.  3a.  which  also  includes 
ray  tracing  synthetics.  The  corresponding  2D  FD  synthetics,  shown  in  Fig.  3b,  are  very 
similar  to  those  generated  by  ray  tracing  except  for  tiie  source  spectrum.  A  noticeable  fea¬ 
ture  in  Fig.  3b  is  the  “ringing”  appearance  of  the  Rg  wavetram  caused  by  the  velocity  gra¬ 
dient  in  the  uppermost  paits  of  the  crust.  Anyway,  the  striking  featiiie  of  the  Fig.  3 
synthetics  is  the  complete  lack  of  significant  coda  waves  which  besides  have  no  counter¬ 
part  in  comparison  to  observational  data  (Fig.  1  and  2).  In  other  words,  for  signal  frequen¬ 
cies  above  1  Hz  homogeneous  crust/lithosphere  models  are  not  tenable  for  computing 
realistic  seismograms. 

Homogeneous  crust  with  sinusoidal  shaped  Moho 

This  model  is  very  simple  except  that  the  Moho  is  given  a  sinusoidal  form  with  ^  =  8 
km  and  Amp  =  1  km.  The  corresponding  2D  FD  synthetics,  dramatically  different  from 
those  in  Fig.  3b,  are  shown  in  Fig.  4.  We  have  here  used  both  P-  and  S-sources  at  a  depth 
of  10  km.  Tlie  striking  feature  of  these  synthetics  is  the  abundance  of  distinct  secondary 
phases  of  types  pP,  SmS,  etc.,  that  is,  crustal  reverberations.  However,  m  contrast  to  coda 
waves  in  real  records,  the  above  secondary  phases  are  highly  correlated  across  line  arrays 
of  apertures  5-10  km.  Naturally,  a  flat  Moho  also  gives  rise  to  correlated  crustal  reverbera¬ 
tions.  Another  interesting  feature,  specific  for  a  corrugated  Moho.  is  P-to-S  converted 
wavelets  (Fig.  4a)  appearing  between  the  major  P-  and  S-arrivals.  In  general,  S-to-P  con¬ 
versions  seem  to  be  more  efficient  than  P-to-S  conversions.  The  Fig.  4a,b  synthetics  also 
demonstrate  that  there  could  be  a  plethora  of  secondary  P-arrivals  and  to  single  out  one  of 
these  as  a  pP-pha.se  for  focal  depth  estimation  is  considered  a  dubious  undertaking. 
Although  using  a  regular,  sinusoidal  .shaped  Moho  interface  iepie.sents  an  oversimplilica- 
tion.  It  serves  to  illustrate  that  even  minor  undulations  of  a  major  structural  discontinuity 
could  profoundly  affect  seismic  recordings  at  local  distances  in  a  nontrivial  way. 
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Random  media  synthetics 


In  this  model  the  von  Kamian  medium  considered  is  of  the  exponential  type  (v  =  0.5) 
with  a  =  0.004  km  and  RMS  =  4  per  cent.  The  a  parameter  value  here  reflects  a  model 
specification  error,  should  have  been  4  km.  but  still  reasonably  realistic  synthetics  were 
produced.  For  the  P-source  synthetics  (Fig.  5a),  the  coda  excitation  is  moderate  relative  to 
the  preceding  P-arrival.  On  the  other  hand,  for  the  first  4  sec  into  the  coda  the  coherency  is 
high,  and  then  drops  sharply  to  a  level  typical  of  observational  data.  The  S-source  synthet¬ 
ics  (Fig.  5b)  appear  somewhat  different  and  on  two  accounts.  Firstly  since  P-excitation  is 
weak,  the  P-phase  amplitude  and  the  coda  level  are  roughly  similar.  Secondly,  the  coda 
coherency  is  rather  high  due  to  a  specific  phase  arrival  at  42  sec  on  the  local  time  scale. 
With  no  ■‘deterministic”  explanation  for  this  pha.se  arrival,  it  is  considered  due  to  fortu¬ 
itous  interferences  of  scattering  wavelets.  Ilic  essence  of  this  experiment  is  that  a  medium 
with  RMS  velocity  perturbations  of  ca  4  per  cent  suffice  for  generating  a  reasonably  strong 
coda.  Observational  studies  also  imply  velocity  perturbations  of  this  order  (4  per  cent)  for 
conelation  distances  in  10-30  km  range  (e.g.,  see  Aki,  1973;  Berteussen  et  al,  1975;  Flatte 
and  Wu,  1988). 

Generalized  random  media  synthetics 

In  this  experiment,  we  use  a  crust/lithosphere  model  which  we  consider  to  be  fairly 
representative  of  the  real  earth.  Firstly,  a  top  crust  low-velocity  layer  (LVL)  of  1 .4  km 
thickness  is  included,  as  this  appears  to  be  a  widespread  continental  structural  feature 
(Ruud  et  al,  1992).  The  medium  randomization  (detailed  in  Table  1)  also  includes  pertur¬ 
bations  of  the  LVL  bottom  interface  and  the  Moho.  'ITie  RMS  velocity  variations  are  4  and 
2  per  cent  in  the  crust  and  below  Moho,  respectively.  Note  that  both  Flatte  and  Wu  (1988) 
and  Charette  (1991),  on  the  basis  of  NORSAR  and  NOR  ESS  scattering  observations, 
favor  relatively  smaller  velocity  perturbations  below  Moho  than  above  Moho  (ca  3  and  2 
per  cent,  respectively).  Anyway,  the  synthetics  for  this  medium  are  displayed  Fig.  6, 
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and  the  following  comments  apply:  Firstly,  the  S-synthetics  in  Fig.  6b  have  some  resem¬ 
blance  to  the  BlisjO  recordings,  although  the  Pn,  Pg  and  PmP  are  somewhat  weak.  The  P- 
source  synthetics  <Fig.  6a)  exhibit  some  resemblance  to  the  profiling  records  displayed  in 
Fig.  2.  The  first-arriving  Pn  is  weak,  followed  by  strong  Pg  and  PmP  phases,  and  then  a 
prominent  and  persistent  coda  wavetrain.  In  an  attempt  to  identify  specific  arrivals  here,  a 
semblance  analysis  was  undertaken.  For  the  P-source  records  (6a),  wavelets  arriving  after 
the  50  sec  time  mark  (semblance  values  0.50  or  greater)  exclusively  exhibited  velocities  in 
the  range  3.0  to  4.0  kms  *.  For  the  S-source  synthetics  (Fig.  6b),  a  semblance  analysis 
gave  similar  velocity  results,  although  the  shear  arrivals  started  at  53  sec  and  semblance 
values  reached  0.8  units.  In  other  words,  Sn-phases  must  be  very  weak  and  the  same 
applies  to  SmSSmSSmS-type  of  phases  with  expected  phase  velocities  in  excess  of  4.4 
kms''.  Seemingly,  many  of  the  late-arriving  wavelets  must  be  caused  by  asymmetric 
reflections  caused  by  the  non-plane  boundaries  of  the  crustal  waveguide. 

As  noted  above,  there  are  several  common  features  between  the  observational  records 
m  Figs.  1  and  2  and  our  synthetics  as  displayed  in  Fig.  6.  However,  there  is  also  a  signifi¬ 
cant  discrepancy  between  these  two  record  sets,  namely,  that  the  synthetic  coda  waves  are 
too  highly  correlated.  The  obvious  explanation  here  is  that  for  computational  reasons  we 
are  restricted  to  2D  modelling  and  a  flat  surface,  thus  ignoring  “out  of  plane”  interfenng 
scattering  wavelets  and  also  scattering  by  topography.  Observationally,  it  is  fairly  easy  to 
demonstrate  the  importance  of  such  contributions,  say  by  array  analysis  of  coda  waves 
(Bannister  et  al,  1991).  Another  problem  is  that  in  the  synthetic  records,  the  Pn-phase 
appears  too  weak  in  comparison  to  Pg-  and/or  PmP-phases.  Tentatively  we  attribute  this  to 
a  lack  of  even  a  weak  sub-Moho  velocity  gradient  -  the  validity  of  this  hypothesis  will  be 
tested  by  computing  synthetics  for  models  incorporating  such  features. 
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6  Discussion 

This  study  was  aimed  at  advancing  our  understanding  of  high-frequency  seismic  w  ave 
propagation  at  local  distances,  thus  exclusively  involving  the  crust/lithosphere  system. 
FTom  numerous  previous  studies  we  have  that  such  media  are  inhomogeneous  and  often 
are  represented  by  random  fields,  where  complexities  can  be  expressed  in  terms  of  a  few 
low  order  statistical  moments  based  on  scattering  parameters  like  correlation  distance  and 
RMS  velocity  perturbations.  Since  traditional  approaches,  say  ranging  from  ray  tracing  to 
perturbation  methods,  are  not  entirely  adequate  for  computing  synthetics  for  inhomoge¬ 
neous  media,  we  have  explored  the  usefulness  of  2D  finite  difference  synthetics.  We 
started  with  simple  structural  models  of  the  kind  used  in  seismic  profiling  studies,  and  not 
surprisingly  found  that  the  corresponding  synthetics  had  little  in  common  with  real  obser¬ 
vations.  Since  Moho  represents  a  first-order  discontinuuy  with  medium  property  contrasts 
of  the  order  of  15-20  per  cent,  perturbing  the  geometry  of  this  interface  as  demonstrated  in 
Fig.  4  strongly  contributes  to  the  number  of  observable  secondary  phases  and  also  the 
coda  per  se.  In  Fig.  5  we  demonstrate  that  an  inhomogeneous  medium  (flat  Moho)  with  an 
RMS  velocity  perturbation  of  4  per  cent  would  generate  a  significant  amount  of  coda 
waves.  Our  most  complex  crust/lithosphere  model  features  both  a  top  crust  low-velocity 
layer,  a  perturbed  Moho  geometry  and  otherwi.se  randomized  velocity  fluctuations  “pro¬ 
duce”  synthetic  seismograms  with  complexities  comparable  to  real  recordings,  although 
the  coda  wave  correlation  is  a  bit  high. 

At  this  stage  of  development,  we  have  not  specifically  aimed  at  quantifying  the  rela¬ 
tive  importance  of  various  scattering  sources.  Even  on  a  very  powerful  IBM  3090  machine 
(IBM  Scientific  Centre  Bergen),  it  takes  days  to  compute  synthetics  like  those  displayed  in 
Fig.  6.  In  the  context  of  event  classification,  geometrical  aspects  of  such  problems  must  be 
explored.  For  example,  to  what  extent  are  generations  of  secondary  phases  and  coda 
waves  dependent  on  focal  depth.  Likewise,  could  Moho  undulations  act  as  a  blocking 
mechanism  for  efficient  crust/lithosphere  wave  propagations?  Such  barriers  are  often 
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attributed  to  mountain  roots,  sedimentary  basins  and  even  ancient  plate  boundaries.  As 
demonstrated  above,  2D  FD  synthetics  are  well-suited  for  providing  an  in-depth  physical 
understanding  of  such  problems. 


7  Concluding  remarks 

In  this  report  we  have  demonstrated  that  2D  FD  synthetics  computed  for  even  simple 
(inhomogeneous)  von  Karman  media  can  produce  quite  complex  seismograms  including 
prolonged  coda  waves.  Although  we  have  not  specifically  attempted  to  evaluate  which  of 
the  medium  parameters  are  dominant  for  generating  coda  waves,  RMS  fluctuations  of  the 
order  of  4  per  cent  and  minor  Moho  perturbanons  are  clearly  important  here.  Also,  a  pure 
P-source  m  a  cracked  medium  would  produce  relatively  strong  Lg-waves  as  typical  of 
many  mining  explosions.  On  the  other  hand,  our  S-source  generates  some  P-waves,  which 
implies  that  S-to-P  conversions  are  more  efficient  than  P-to-S  conversions.  In  all  our  syn¬ 
thetic  records,  the  Pn-phase  is  relatively  weak,  which  we  take  to  imply  that  sub-Moho 
constant  medium  velocity  probably  is  not  a  good  approximation  of  the  real  earth. 

Finally,  the  advantage  of  FD  synthetic  seismogram  analysis  is  an  improved  physical 
insight  in  complex  media  like  the  crust/lithosphere  system.  This  in  lurn  is  of  importance 
for  the  design  of  event  classification  criteria  simply  because  the  paili  effect  and  focal  depth 
would  be  dominant  factors  at  local/regional  distance  ranges. 
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Table  1 :  Crust/lithosphere  model  used  for  computing  the  synthetics  for  P-  and  S -sources  at 
a  depth  of  12  km  (waveforms  are  displayed  in  Fig.  6).  H  is  layer  thickness.  S-veloc- 
ities  were  calculated  from  the  P-velocities  (Pvel)  given  a  Poisson  ratio  of  0.25.  For 
the  von  Karman  random  media,  v  is  the  order  and  a  is  the  correlation  distance  in  km. 
The  RMS  fluctuation  is  given  in  km  for  the  corrugated  interfaces  and  in  per  cent  for 
velocities. 


112 


1985/179/15;42;40.0  NORESS  AO  SP  Azi  240  Dist  300  km 


113 


Fig.  1:  NORESS  recording  of  a  BiSsjo-explosion.  The  outstanding  scattering-related  foa- 


REDUCED  TIME  (sec)  (8  kms‘ 


Fip.  3  Comparison  between  ray  uacinp  and  2D  Fd  synthetics  for  a  “High  Velocity  Lower 
Cnist"  niixlcl  taken  from  Ilyndman  and  Klcm(reier  (198V),  For  such  a  simple  homo 
geneous  model,  the  diffeiences  between  the  two  .sets  of  synthetic  iccoids  aie  small. 


A  peculiar  feature  in  the  I  D  -s'  li.eiics  (I'ig.  h)  is  the  app.iient  ringing  in  the 
Rayleigh  waves  apparently  c.iu-  cd  by  the  velocity  giadienis  in  the  upix'r  ciaist 
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Fig.  4,  2D  FD  .synthelics  for  a  simple  crustal  model,  but  with  a  conugated  Moho  (X  =  8 
km,  amp  =  1  km).  Both  P-  and  S-souices  are  used,  and  the  focal  depth  was  lO  km. 
The  lower  parts  of  figure  a)  and  b)  arc  semblance  analysis  results  of  the  respective  P- 
and  S-souri  e  synthelics.  .Some  of  the  secondary  anivals  arc  likely  to  have  a  scatter¬ 
ing  origin  (poor  semblance)  and  al.so  that  S-to-P  conversions  appear  to  be  more  effi¬ 
cient  than  P-to-S  conversions. 
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Fig.  6:  P-  and  S-source  synthetics  for  the  medium  detailed  in  Table  1 .  Source  depth  was  1 2 
km  and  epicentral  distance  was  230  km.  Below  the  respective  P-  and  S-source  traces 
the  corresponding  semblance  analysis  results  are  given.  The  outstanding  feature 
here  is  that  all  prominent  phase  arrivals  after  48  sec  (P)  and  52  sec  (S)  have  phase 
velocities  in  the  range  3-4  kms  ’. 
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ABSTRACT 

The  problem  of  discriminating  between  earthquakes  and  underground  nuclear  explosions  is 
formulated  as  an  exercise  in  pattern  recognition  approach  analysis.  An  advantage  of  our  procedure 
is  flexibility,  by  combing  both  adaptive  noise  suppression  and  event  classification  incorporating 
feature  selection  criteria. 

The  analysis  has  been  applied  to  a  learning  set  of  44  nuclear  explosions  (8  test  sites)  and  35 
earthquakes  in  Eurasia  recorded  at  the  NORESS  array  (Fig.  1).  The  signal  features  considered  were 
the  normalized  power  in  8  spectral  bands  in  the  0.2-5.0  Hz  range  of  the  P-wave  (6  sec)  and  the  P- 
coda  (30  sec).  Physically,  it  means  that  we  exploit  potential  differences  in  the  shape  of  earthquake 
and  explosion  spectra,  respectively.  Other  features  included  are  peak  P  and  P-coda  amplitude 
frequencies  and  relative  P/P-coda  power.  These  19  features  were  extracted  either  from  conventional 
array  beam  traces  or  the  optimum  group  f  iltered  traces  (OGF-removal  of  coherent  low-frequency 
noise).  Using  the  feature  selection  algorithm,  based  on  estimates  of  the  expected  probability  of 
misclassification  (EPMC),  only  2  to  4  features  were  needed  for  optimum  discrimination 
performance.  The  dominant  features  were  coda  excitation  and  P-  and  P-coda  power  at  lower  signal 
frequencies.  Furthermore,  feature  parameters  extracted  from  the  OGF  traces  had  a  slightly  better 
performance  in  comparison  to  those  extracted  from  beam  traces.  Finally,  there  were  no 
misclassifications  for  OGF-derived  features  when  the  explosion  population  was  limited  to  E. 
Kazakh  events,  while  including  events  from  the  other  test  sites  lead  to  a  decrease  in  discrimination 
power. 
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INTRODUCTION 


The  problem  of  distinguishing  underground  nuclear  explosions  from  natural  earthquakes  using 
seismic  data  has  been  studied  (or  a  long  time.  In  fact,  it  dates  back  to  19  Sep  1957  as  on  that  day 
a  nuclear  explosion,  code-named  Rainier,  was  detonated  under  the  Nevada  desert.  The  principal 
goal  of  the  experiment  was  to  explore  the  ability  of  an  underground  test,  unhampered  by  weather 
and  concerns  over  radioactive  fallout,  to  fulfill  all  the  needs  of  a  nuclear  test  program.  In  this 
respect  Rainier  was  a  great  success.  The  resulting  seismological  data  were  studied  intensely  in 
scientific  and  political  circles,  setting  a  pattern  that  still  prevails.  The  introduction  of  seismology  to 
the  international  political  arena  took  place  a  year  later  when  scientific  experts  from  the  UK,  USA 
and  USSR  met  in  Geneva,  in  Aug  1958,  to  design  a  seismic  verification  system  as  part  of  a 
comprehensive  Nuclear  Test  Ban  Treaty  (CTB).  Scientific  experts  (now  from  some  26  countries) 
still  meet  regularly  in  Geneva  to  di.scuss  the  design  of  such  a  system  for  a  potential  CTB.  Source 
identification  (SI)  by  seismological  means  remains  problematic  for  small  events. 

In  order  to  avoid  the  once  troublesome  issue  of  in-country  operation  of  non-national  seismic 
stations,  the  SI  research  in  the  1960  and  1970-ties  focused  on  observations  in  the  tcicseismic 
window.  The  most  successful  criteria  for  seismic  source  identification  were  spectral  ratio  variants, 
mainly  between  non-overlapping  frequency  bands  for  the  P-signal  itself,  or  the  relative  signal 
excitation  at  1  sec  (P-wave)  and  20  sec  (surface  waves).  The  latter  is  often  denoted  the  m^iM, 
(body  wave  versus  surface  wave  magnitudes)  discriminant.  Another  commonly  used  discriminant 
was  the  so-called  complexity  tied  to  the  ratio  of  P  coda  RMS  in  two  consecutive  windows  of 
lengths  around  5  sec  and  30  sec  respectively  (Dahlman  and  Israelson,  1977;  Douglas,  1981).  A 
variant  of  the  complexity  often  denoted  the  coda  discriminant,  was  introduced  by  Tjdsthcim  (1975, 
1978)  and  Sandvin  and  Tjdstheim  (1978)  using  autoregressive  spectral  coefficients  as  signal 
attributes  in  combination  with  more  advanced  discrimination  statistics.  A  vast  literature  exists  on 
the  seismic  source  identification  problem  and  the  theoretical  foundation  for  the  mentioned 
discriminants;  for  example,  sec  Bolt  (1976),  Dahlman  and  Israelson  (1977),  Husebye  and 
Mykkcltvcit  (1981),  Kerr  (1985).  Prc.ss  (198.5),  Tsipis  et  al  (1986)  and  the  many  articles  contained 
in  these  books.  Other  references  arc  Evemden  (1977);  Evemden  and  Kohler  (1979),  and  Evemden 
et  al  (1986),  Blandford  (1982),  and  Pomeroy  ct  al  (1982). 
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The  above  discriminants,  and  in  particular  the  m^:  M,  one,  arc  ctficient  in  terms  of  few  failures, 
but  often  fail  for  events  with  m,,  magnitudes  below  4.5  -  5.0  units  due  to  detoriating  signal -to-noise 
(SNR)  ratios.  Thi.s,  seemingly  in  combination  with  a  lack  of  novel  ideas  for  source  discriminants 
in  the  tclescismic  window,  is  taken  to  explain  the  relatively  few  discrimination  papers  published 
during  the  1980s.  On  the  other  hand,  during  this  decade  seismic  instrumentation  and  array  design 
improved  considerably.  It  suffices  here  to  mention  the  advent  of  the  24  bits  digitizer,  increased 
bandwidths  (20  Hz  or  higher)  and  the  introduction  of  new  arrays  like  NORESS  for  event 
monitoring  at  local  distances.  In  these  ranges  the  event  detectability  is  likely  to  be  around  2.5  to 
3.0  magnitude  units  at  the  90  per  cent  level  (Sereno  et  al,  1991)  and  the  corresponding  SI 
performance  is,  by  rule-of-thumb,  likely  to  be  in  the  3.0  to  3.5  magnitude  range.  Obviously,  a 
signal  must  be  relatively  stronger  in  ease  of  classification  as  compared  to  detection  analysis.  The 
high  quality  data  now  at  hand  from  modem  arrays  and  3-componcnt  stations  (also  USSR 
deployment)  have  again  triggered  interest  in  seismic  discrimination  problems  (Baumgardt  and 
Young,  1990).  Also  new  concepts  have  been  introduced,  namely,  the  so-called  artificial  or  trained 
neural  networks  technique  (e.g.  see  Dowla  et  al,  1990;  Dysart  and  Pulli,  1990). 

In  this  study  we  will  address  the  problem  of  tcleseismic  source  discrimination  and  explore  the 
potential  of  the  spectral  ratio  and  complexity  discriminants.  This  may  seemingly  be  a  step 
backward,  but  it  is  easy  to  argue  tiiat  it  is  not  so.  We  use  data  from  the  NORESS  array,  which  has 
an  excellent  detectability  for  events  in  parts  of  Eurasia  (Ringdal,  1990),  and  besides  many  regions 
arc  still  lacking  adequate  seismograph  network  coverage.  Also,  use  of  coda  waves  may  prove 
instructive  for  similar  di.scriminants  at  local  and  regional  distances  which  presently  arc  tied  to  the 
spectral  content  of  crustal  phases  like  Pg,  Pn,  Sn,  Sg  (Eg).  We  add  here  that  the  recent  Taylor  and 
Marshall  (1991)  discrimination  study  was  based  on  UK-type  array  recordings  (YKA,  EKA,  GBA, 
WRA)  of  Eurasian  events  also  in  the  teleseismic  window. 


EVENT  SELECTION  AND  NORESS  RECORD  PREPROCESSING 

The  presumed  earthquakes  (PDE  and  ISC  listings)  and  presumed  underground  nuclear  explosions 
(NO.  SAR  listings)  used  in  this  study  are  given  in  Tables  1  and  2.  According  to  P.  Richards  (pers. 
comm.)  such  event  listings  are  not  always  foolproof,  hence  we  have  used  the  word  "presumed" 
here.  The  latter  comprises  alt  Soviet  explosions  recorded  at  the  NORESS  array  which  became 
operational  in  late  fall  1984  except  for  5  presumed  explosions  for  which  no  data  were  available  to 
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US  (details  in  Tabic  1  caption).  Likewise,  data  were  lacking  for  14  presumed  earthquakes  for  which 
the  epicenters  were  within  the  "framed"  area  in  Fig.  1.  The  event  epicenters  are  depicted  in  Fig.  1 
as  well.  Most  of  the  explosions  (32  out  of  44)  stem  from  the  E.  Kazakh  (Semipalatinsk)  test  site 
and  hence  the  earthquakes  in  the  m,,  range  4.0  to  6.0  were  confined  to  the  same  general  area.  The 
reason  was  that  the  record  differcncas  .should  mainly  be  attributable  to  source  parameters  and  no. 
reflect  significant  differences  in  propagation  paths.  The  other  test  sites  arc  in  ascismic  areas  for 
which  sufficient  earthquake  recordings  arc  not  available.  However,  the  outlaying  explosions  were 
included  for  a  check  on  discriminant  robustness.  The  "choice"  of  our  explosion  and  earthquake 
populations  reflects  source  identification  outcomes  as  reported  by  the  ISC  or  PDE  (NEIC/USGS) 
agencies.  There  is  hardly  any  other  way  for  "event"  selections;  in  case  of  gross  errors  here  we  aim 
to  identify  possible  spurious  events.  Also,  the  explosions  are  presumed  to  be  nuclear  ones  for  the 
simple  reason  that  besides  being  located  in  well-established  test  site  areas,  chemical  explosions  on 
land  are  hardly  ever  recorded  at  teleseismic  ranges.  The  small  NORESS  array  of  aperture  3  km  has 
an  excellent  event  detectability  for  parts  of  Eurasia.  With  a  40  Fiz  sampling  the  array’s  bandwidth 
is  20  Hz.  A  simple  and  efficient  scheme  for  SNR  enhancement  is  delay-and-sum  processing  or 
beamforming  which  is  most  efficient  in  the  2-8  Hz  band  (e.g.  see  Birtill  and  Whiteway,  1965; 
Ingaie  et  al,  1985;  and  Husebyc  and  Ruud,  1989).  At  lower  frequencies  (below  say  2  Hz)  the 
wavelengths  of  microseisms  are  of  the  same  order  as  the  array  aperture,  and  hence  strong 
correlation  in  the  noise  across  the  array  is  often  observed,  in  such  cases  maximum  likelihood 
schemes  are  very  efficient  in  suppressing  correlated  noise,  as  demonstrated  by  Ingate  et  al  (1985). 
Recently,  even  more  advanced  methods  have  been  introduced  by  Kushnir  et  al  (1990),  which  are 
extensively  used  here  for  suppressing  low  frequency  propagating  noise.  The  efficiency  of  the  this 
scheme,  commonly  denoled  the  Optimal  Group  Filtering  (OGF)  technique  is  demonstrated  in  Fig.  2 
and  also  in  Fig.  9.  Naturally,  removal  of  low  frequency  noise  is  important  as  both  theoretical  and 
observational  studies  demonstrate  that  part  of  the  discrimination  power  is  vested  in  the  low/high 
frequency  bands  of  the  P-signal  (Evemden  et  al,  1986;  Taylor  and  Marshall,  1991).  In  the 
Tjdstheim  and  Husebye  (1976)  and  Sandvin  and  Tjostheim  (1978)  studies,  the  noise  suppression 
was  by  beamforming,  which  for  the  large-aperture  NORSAR  array  was  efficient  due  to  the  much 
larger  sensor  spacing  here.  In  addition,  noise  spectral  estimates  based  on  the  preceding  noise  were 
subtracted  from  the  P-signal  spectral  estimates  but  with  marginal  improvements.  We  also  tried  this 
kind  of  noise  subtraction  and  with  a  similar  outcome. 

Evemden  (1977)  advocates  the  use  of  P-signal  frequencies  up  to  9  Hz  in  teleseismic  discrimination 
studies,  while  we  restricted  the  analysis  to  5  Hz.  Our  rationale  being  that  there  is  not  much  signal 
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energy  above  5  Hz  and  besides  cultural  sources  like  local  quarry  blasting,  mining  explosions  and 
fast  running  machinery  (saw  mills)  could  easily  bias  the  observational  data.  Also,  our  discriminant 
parameters  were  extracted  both  from  single  P-beam  traces  and  from  OGF  traces  in  order  to  have 
observations  for  judging  the  relative  importance  of  the  latter. 


SIGNAL  ATTRIBUTES  -  CLASS  OF  WAVEFIELD  PARAMETERS  FOR  DISCRIMINATION 

As  mentioned,  the  most  powerful  discrimination  prameters  are  related  to  differential  signal 
excitation  in  different  frequency  bands  for  earthquakes  and  explosions  respectively.  P-wave 
prameters  are  an  obvious  choice  here,  because  this  phase  is  most  easily  detected.  P-coda  waves 
are  intesting  as  they  not  only  reflect  the  source  but  also  the  source  location  within  the  crust/uppr 
mantle.  For  example,  coda  excitation  and  duration  are  far  less  efficient  for  surface  explosions  and 
deep  earthquakes  relative  to  shallow  and  intermediate  depth  earthquakes  (Dainty,  1990).  Rayleigh 
waves  are  not  considered  simply  because  they  become  embedded  in  background  noise  for  event  m^ 
magnitudes  at  4.5  -  5.0  or  below.  Another  disadvantage  is  that  of  extensive  interference  of  surface 
waves  from  other  events  are  likely  to  occur  (Marshall  and  Douglas,  1985).  The  essence  of  the 
above  discussion  is  that  our  class  of  discrimination  prameteis  is  tied  to  the  P-signal  and  its  coda 
as  illustrated  in  Figs.  2,  3  and  4.  The  choice  of  window  lengths  of  6  sec  and  30  sec  reflects  that 
most  of  the  desirable  information  from  teleseismic,  but  also  regional  events  are  accumulated  in  the 
first  3-6  sec  of  the  P-wave  and  then  the  subsequent  10  -  30  sec  coda  waves  where  scattering 
contributions  arc  still  significant  (e.g.,  sec  Tjdstheim,  1981;  Dainty  and  ToksCz,  1990;  Bannister  et 
al,  1990). 

Processing  details  were  as  follows.  For  each  event  5  minutes  of  recordings  for  all  25  vertical 
NORESS  sensors  were  extracted.  The  first  2  minutes  of  pure  noise  were  used  for  estimating  OGF- 
filter  coefficients  for  the  slowness  and  azimuth  of  the  individual  events.  After  the  filtering  was 
prformed,  amplitude  spectra  (FFT)  were  calculated  for  the  P-signal  (6  sec)  and  the  P  coda  (30 
sec).  The  power  spetra  for  the  non-overlapping  8  bands  specified  in  Figs.  2  and  3  were  obtained 
by  simple  averaging  of  spctral  squared  amplitudes.  Since  the  events  used  in  analysis  have  widely 
different  magnitudes,  all  power  spetra  were  normalized  by  their  maximums.  Physically,  this  means 
that  ptential  spctral  shap  differences  are  exploited  for  event  discrimination  purpses.  The  final 
feature  prameter  values  were  obtained  by  taking  the  logarithm  of  the  normalized  spectral  values. 
As  shown  in  Fig.  3  we  have  8  feature  prameteis  for  both  the  P  signal  and  the  P  coda  for  the  same 
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set  of  frequency  bands.  Additional  parameters  introduced  (nos.  17  and  18)  are  peak  spectral 
frequencies  and  finally  the  19th  parameter  being  the  ratio  of  P/P  coda  spectral  maxima.  Note  that 
the  last  parameter  is  close  to  the  classical  complexity  definition.  Also,  from  past  studies  the  largest 
spectral  differences  between  Eurasian  earthquakes  and  explosions  appear  to  be  in  the  2.6-3.2  Hz 
band  according  to  Dahiman  and  Israclsson  (1977),  and  0.5-1.0  Hz  and  2.0-3.0  Hz  bands  according 
to  Taylor  and  Marshall  (1991).  Finally,  normalized,  average  power  spectra  for  the  8  frequency 
bands  detailed  in  Fig.  3  are  for  all  events  analyzed  shown  in  Fig.  4. 

DISCRIMINATION  APPROACH 

The  problem  of  discriminating  between  earthquakes  (EQ)  and  underground  nuclear  explosions 
(NE)  is  of  a  general  nature.  Hence,  there  are  many,  often  complex  classification  algorithms  which 
are  of  potential  relevance  in  this  context.  Here,  we  may  differentiate  between  non-statistical 
approaches  like  those  of  "neural  networks"  (e.g.,  see  Dowla  et  al,  1990)  and  the  "Koia-algorithms" 
of  Gelfand  et  al  (1976),  which  contrast  with  a  statistical  approach  to  pattern  recognition 
(Tjdstheim,  1981).  Also,  the  seismic  discrimination  undertaking  is  a  two-stage  process:  firstly, 
relevant  discrimination  parameters  must  be  defined  and  extracted  from  the  records,  and  secondly, 
descision  rules  (discriminators)  must  be  introduced  to  ensure  proper  event  classification. 

At  the  outset  of  this  study,  we  discussed  rather  extensively  among  ourselves  which  discrimination 
approach  would  be  best  for  cla.ssification  of  Eurasian  events  as  recorded  by  NORESS.  The 
importance  of  (he  physical  aspect  of  the  problem  at  hand  was  duly  recognized,  that  is,  the 
extracted  discrimination  parameters  must  be  seismologically  relevant.  For  example,  Tjdstheim  and 
Sandvin  (1978)  found  that  P-wave  and  coda  autoregressive  spectral  parameters  were  efficient 
discriminants,  although  their  scismological  relevance  is  not  obvious.  Furthermore,  the  persistence 
of  the  seismic  source  identification  problem  is  that  a  unique  set  of  discriminants  for  weak  events 
apparently  does  not  exist,  so  we  wanted  to  have  the  ability  to  decide  statistically  which  ones  of 
such  parameters  would  be  most  informative  for  a  given  area  or  region.  It  is  here  natural  to  use  the 
probability  of  misclassification  (PMC)  as  a  measure  for  discrimination  power  and  on  this  basis  to 
select  most  significant  feature  parameters.  In  short,  we  decided  to  use  a  statistical  pattern 
recognition  approach  which  is  well  suited  for  PMC  estimation  when  the  number  of  observations  is 
small  and  the  number  of  classification  parameters  is  relatively  large. 


The  recognition  problem  in  our  case  may  be  formulated  as  follows.  The  data  set  consists  of  three 
sets.  The  first  two  are  the  sets  R,  of  explosions  (NE)  spectral  parameters  vectors  ri^‘\  and  Rj 
earthquake  (EQ)  spectral  parameters  vectors  rf^: 
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number  of  observations  in  the  j-th  learning  set  (earthquakes  or  explosiciis),  j=l,2 
values  of  classification  features  (spectral  parameters  P  and  P-coda),  k=l,n 
number  of  classification  features. 


R,  and  Rj  form  the  learning  material  for  the  NE  and  EQ  classification  features.  The  third  set  is 
the  vector  X  containing  discrimination  parameters  for  an  "unknown"  event,  i.e.,  the  vector  being 
classified  with  no  prior  knowledge  as  to  its  source  identification.  On  the  basis  of  the  learning 
samples  R,  and  Rj,  we  must  make  a  decision  on  whether  X  belongs  to  the  first  or  the  second  class. 
This  is  done  by  using  a  decision  (discriminator)  function  g(X,  R„  R^).  The  equivalent  Bayesian 
rules  are:  If  g(X,  R,,  R2)  <  C,  the  hypothesis  HI  is  correct  and  X  belongs  to  the  Orst  class;  if  g(X, 
Ri,  Rd  >  C,  the  hypothesis  H2  is  correct  and  X  belongs  to  the  second  class.  Here  C  is  a  constant, 
that  is,  a  threshold  value.  Using  this  decision  rule,  there  are  errors  of  two  kinds:  the  vector  X 
belongs  to  the  second  class  (the  hypothesis  H2  is  correct),  while  it  is  assigned  to  the  Gist  class. 
Errors  of  the  second  kind  represent  the  reverse  situation.  Errors  occur  because  the  observations  R„ 
Rj  and  X  are  random,  measurement  errors  and  the  random  nature  of  the  discrimination  features 
themselves.  The  classification  errors  of  the  first  and  second  kind  are  described  by  the  probabilities 
PI  and  P2.  It  is  also  assumed  that  a  priori  probabilities  q^;  j=l,2  of  the  vector  X  belonging  to 
the  j-th  class  are  known,  and  the  total  error  probability  P„  =  qi  PI  <iJ?2  can  then  be  evaluated. 
In  our  case  the  sizes  of  the  learning  samples  (35  EQ  and  44  NE)  are  comparable  to  the  dimension 
n  of  the  features  space  (19  spectral  parameters).  In  such  cases,  it  is  difDcult  to  resolve  the  basic 
classiGcation  problem  ~  to  select  the  optimum  decision  rule  and  evaluate  its  corresponding  error 
probabilities.  Since  the  distribution  patterns  corresponding  to  the  Grst  and  second  class  are 
unknown,  it  is  theoretically  impossible  to  construct  a  uniformly  optimum  decision  rule,  which  in 
all  cases  will  yield  the  least  probability  P„  of  missclassiOcation  (Kushnir  et  al,  1986).  For  this 
reason  a  large  number  of  different  rules  with  different  "good"  properties  is  used  in  practice.  The 
criterion  for  choosing  a  particular  rule  is  the  requirement  of  minimizing  the  PMC.  Qther  aspects, 
such  as  computational  efOciency  and  the  possibility  of  analytical  evaluation  of  the  PMC  associated 
with  a  particular  rule  should  also  be  considered. 
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From  this  point  of  view,  when  the  R,  and  Rj  distributions  arc  norma',  the  linear  discriminant 
function  (LDF)  is  optimal  in  some  statistical  sense  (Anderson,  1958). 

L(X,R„Rj)  =  (X-0.5(r,  +  r^))*  S'Vi  -  r^)  <  C  ,  (2) 

where 

Ni 

r^  =  1/N  Z  r,«  ;  j=l,2 
i=l 

S  =  I(l/(Nj  +  Nj -1)]  (A,  +  Aj) 
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A^=l/(Nj-l)  E  (r.®  -  r^r.®  -  r/ 
i=l 

X  -  vector  being  classified, 

S  -  covariance  matrix  of  learning  sets, 

Tj  •  average  vector  of  j-th  learning  set, 

C  -  threshold;  above  is  EQ,  below  is  NE 

Nj  -  number  of  observations  in  the  j-th  learning  set 

•  -  matrix  transpose. 

Note,  the  LDF  expression  in  eq.  (2)  is  often  denoted  the  Fischer  discriminant  (Anderson,  1958). 
When  the  covariance  matrices  of  learning  samples  are  unequal  and  the  distributions  are  norma)  the 
quadratic  discriminator  (QDF)  may  be  used: 

Q(X,  R„  Rj)  =  (X  -  r,)*  S,  '(X  -  r,)  -  (X  -  r^)*  Sj'(X  -  r^  -ln(q,  /  q^)  >  C  ,  (3) 

where 

Sj  -  covariance  matrix  of  j-th  learning  set,  j=l,2, 
qj  -  a  priori  probability  of  j-th  learning  set. 


I£rror  probability  estimation 


In  the  ease  of  Tixed  learning  samples  R,  and  R.  the  decision  rule  g(X.  R„  R,)  is  a  function  of 
only  one  vector  -  the  vector  l)cing  classilied  X.  I  he  PMC^i  wcurring  upon  llxing  R,  and  Rj  arc 
called  conditional  probabilities.  I  hcse  piobabilities  determine  the  power  of  the  rule  "trained"  once 
and  then  used  repeatedly  alterwards.  Conditional  piobabilities  are  random  and  on  changing  the 
learning  samples  they  can,  in  principle,  vary  within  wide  limits.  The  mean  values  of  conditional 
probabilities  with  respect  to  the  distribution  of  the  learning  samples  arc  called  the  expected 
probabilities  of  misclassification  (EPMC).  Estimating  the  EPMC  which  takes  into  account  both  the 
random  nature  of  the  learning  samples  R,  and  R,  and  the  vector  examined  X  is  problematic,  even 
under  the  simple  assumptions  that  the  classes  have  normal  distributions.  For  LDF  an  exact 
expression  of  EPMC  is  too  complex  for  practical  purposes.  However,  an  approximate  expression  of 
the  EPMC  provides  sufficient  accuracy  for  the  parameter  values  used  in  practice.  The  principal 
method  of  deriving  these  approximate  exprc.ssions  for  error  probabilities  is  that  of  assymptolic 
representations  of  the  distributions  of  di.scriminator  statistics.  Of  major  interest  in  seismic 
discrimination  analysis  is  the  case  when  N  (number  of  events  in  learning  sets)  and  n  (number  of 
features)  are  of  the  same  order  and  have  double  digit  values.  For  LDF  (Linear  Discriminant 
Function)  the  asymptotic  equation  suggested  by  Kolmogorov,  where  the  enor  probabilities  are 
investigated  when  both  n  and  N  approach  infinity  and  n/N  is  constant,  remains  approximately 
valid.  Deev  (1970)  has  shown  that  in  Kolmogorov’s  asymptotic,  the  distribution  of  the  LDF 
linear  discriminator  (L  in  cq,  (2))  is  asymptotiailly  normal  with  mean  M  and  variance  V: 

Mj  =  (2(N-l)/(2N-l-n)|  |(-l)'  D’/2|.  j=1.2 

=  |(N-1)(2N-I)(2N+I)  /  (2N-n-l)(2N-n)n|  |D'  +  2n/N),  (4) 

if  N  =  N,  =  Nj, 

where  D^  =  (r,  -  rj'  S  '(n  -  t’)- 

From  (4),  the  following  approximate  expression  for  the  EPMC  of  the  LDF  can  be  written: 

PI  =P{gsC|M2}=G((C-M,)/V) 

P2  =  P{gsC|Hl}=G((Mj-C)/V),  (5) 
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where 
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G(y)  =  Ih/ht  J  c  dx 
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Raudis  and  Pikyalis  (1975)  have  shown  that  the  relative  error  in  evaluating  the  probability  of 
misclassification  yielded  by  asymptotic  expression  (5),  as  compared  with  the  error  computed  from 
an  exact  expression,  is  not  greater  then  few  per  cent  in  the  common  observations  ranges  of  N  and 
n.  The  event  classification  strategy  as  adopted  in  this  study  is  visualized  in  Fig.  5. 

The  LDF  is  used  for  cases  when  the  covariance  matrices  of  the  learning  sets  being  classified  are 
the  same.  The  underlying  assumption  here  is  not  always  valid,  so  we  also  considered  the  more 
complex  quadratic  discriminator  (QDF  -  eq.  (3)),  which  is  useful  when  the  R,  and  R^  matrices  are 
different.  An  expression  for  the  EPMC  for  QDF  has  been  obtained  by  Levin  and  Troitskii  (1970), 
and  will  not  be  given  here.  Anyway,  in  terms  of  the  data  used  in  analysis  in  this  study,  we  found 
that  the  linear  discriminator  gave  fewer  misclassification  errors  than  the  QDF.  We  note  that  even 
under  more  adverse  conditions  with  non-equal  covariance  matrices  and  the  sample  size  N 
comparable  to  the  number  of  features  n,  the  LDF  is  still  preferable  to  the  quadratic  discriminator 
(e.g.,  see  Pisarenko  et  al,  1983).  This  is  mainly  due  to  a  greater  number  of  unknown  QDF 
parameters  in  comparison  to  the  corresponding  LDF  parameters  for  small  learning  sets. 

Broad  applications  of  LDF  in  various  classification  problems  have  confirmed  its  efficiency  for 
small  sizes  N1  and  N2  of  learning  samples  and  for  distributions  that  differ  from  normal  (Azen  et 
al,  1975;  Weber  et  al,  1986;  Tsvang  et  al,  1986).  In  essence,  the  advantage  of  LDF  is  its 
simplicity,  and  the  existence  of  a  method  for  calculating  EPMC. 

Selection  of  most  informative  features 

As  mentioned  above,  the  linear  discriminant  function  has  a  good  perfomance  when  the  sizes  of 
learning  sets  are  small.  LDF  can  be  applied  not  only  to  the  total  number  of  discriminant  features 
in  eq.  (1)),  but  also  to  vectors  made  up  of  a  subset  of  the  discriminants.  Various  values  of  the 
expected  probability  of  misclassification  (eq.  (5))  will  be  obtained.  For  a  certain  optimum  subset  of 
features,  the  EPMCs  for  a  small  size  of  the  learning  set  may  be  much  smaller  than  for  the  total 
system  of  features.  This  unexpected  "multidimensional  effect"  is  explained  as  follows;  the  EPMC 
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for  LDF  arc  expressed  l)y  (cqs.  (4)  and  (5)),  which  depends  on  ihe  dimcnsionalily  of  the  space  of 
features  n,  the  sizes  of  learning  sets  N,  and  Nj,  and  the  Mahalonobis  distance  D  (eq.  (4)).  The 
latter  is  a  non-decreasing  function  of  the  number  of  features  n.  In  practice,  this  function  will  for 
increasing  n  always  tends  to  a  limit,  because  the  information  value  of  features  is  exhaustible 
(usually,  there  arc  few  informative  parameters!).  Correspondingly,  the  EPMC  should  pass  through  a 
minimum,  bccau.se  less  and  Ic.ss  information  is  being  added,  while  the  noise  is  being  increased. 
Formally,  the  existence  of  such  a  minimum  follows  from  cxpic.;sions  (4)  and  (5).  The  natural 
strategy  would  be  to  take  only  the  set  of  features  which  rcprc.sents  the  minimum  of  the  EPMC.  In 
order  to  make  use  of  this  idea,  one  must  be  able  to  arrange  the  available  features  in  a  sequence 
such  that  increasing  the  set  of  features  by  one  feature  at  a  time  would  produce  the  maximum  rate 
of  growth  of  the  function  D\n).  One  possible  way  of  ranking  the  features  is  the  following:  at  each 
iteration  step  K  the  optimum  subse'  of  K-1  features  previously  selected  is  increased  by  adding  the 
one  feature  from  among  the  remaining  set,  which  yields  the  largest  difference  D\K)-D\K-1).  The 
features  arc  thus  ranked  in  significance  and  in  relation  to  the  Mahalonobis  distance  D  as  a  function 
of  features  n.  Then  the  EPMC  are  computed  and  plotted  versus  ranked  features.  The  optimal  subset 
of  features  corresponds  to  the  minimum  of  EPMC,  as  illustrated  in  Fig.  6. 

This  optimal  parameters  algorithm  in  combination  with  the  LDF  (linear  discriminator)  exhibits 
good  perfomance  in  many  geophysical  applications  when  the  sample  sizes  are  small  (Weber  et  al, 
1986;  Tsvang,  1986;  Tsvang  et  al,  1986;  Avsjuk  et  al,  1988).  That  is  why  the  LDF  for  optimal 
parameters  subset  was  chosen  in  our  case  (35  EO  and  44  NE).  In  practice,  when  the  distributions 
of  the  classes  (the  learning  set  ob.scrvation.s)  arc  unknown,  the  results  of  discrimination  are  usually 
obtained  using  so-called  non -parametric  methods.  Such  an  approach  is  very  common  in  seismic 
event  di.scrimination  analysis  and  is  illustrated  in  Fig.  7. 

When  the  data  set  is  large  enough,  the  simplest  way  of  checking  reliability  of  a  chosen 
discriminator  is  by  using  part  of  the  set  for  learning  and  the  remainder  for  testing.  Usually  when 
the  sample  size  is  small,  one  uses  for  the  classiFication  test  the  learning  set.  It  is  so  called  the 
Reclassification  method.  Here,  training  of  the  discriminator  is  carried  out  with  respect  to  the  entire 
size  of  the  learning  samples.  Then  all  events  of  these  samples  arc  sequentially  classified.  It 
provides  usually  too  optimistic  discrimination  results  (Pisarenko  et  al,  1981)  because  learning  and 
testing  just  coincide  (.sec  also  Table  3).  To  provide  more  confident  decision  we  used  also  the 
so-called  Jack-knife  method,  which  allows  the  creation  of  an  independent  data  set  for  testing  the 
discriminator  without  decreasing  the  learning  data  set.  Here  learning  and  testing  are  carried  out 


129 


repeatedly:  during  each  learning  sequence  one  element  is  omitted  from  the  entire  sample  and  in 
turn  used  for  testing.  The  total  number  of  errors  is  then  summarized.  In  essence,  all  elements  of 
the  sample  are  used  both  for  learning  and  testing,  but  the  element  being  used  for  recognition  are 
statistically  independent  of  the  learning  set.  Finally,  it  should  be  noted  that  Jack-Knife  method  may 
also  be  used  for  estimating  the  conditional  probability  of  misclassification  (not  the  EPMC,  see 
section  "Error  probability  estimation"). 


RESULTS 

Perhaps  the  most  important  aspect  of  the  seismic  source  identification  problem  is  that  of  finding 
proper  discriminant  parameteis.  In  practice,  this  has  proved  rather  difficult  since  source  locations 
within  the  crust  and  the  propagation  path  to  the  receiver(s)  are  of  significance  in  this  context.  In 
our  case  we  selected  19  discriminant  candidates  (Fig.  3)  and  have  used  the  EPMC  measure  for 
ranking  the  relative  importance  of  these  parameters  as  illustrated  in  Figs.  6  and  7.  In  these  figures 
we  differentiate  between  two  populations,  namely,  ail  events  (44  NE  and  35  EQ)  and  only 
Semipalatinsk  explosions  and  all  earthquakes  (32  NE  -t-  35  EQ)  for  OGF  and  BEAM  traces, 
respectively.  Parameter  no.  19  (Fig.  3)  are  clearly  the  most  dominant  for  all  cases,  and  thus 
demonstrate  the  relevance  of  the  classical  complexity  discriminant.  There  appear  to  be  some 
differences  between  the  OGF  and  BEAM  trace  derived  parameters,  as  in  the  former  case  the 
optimum  performance  (EPMC  minimum)  is  tied  exclusively  to  the  coda  parameters.  For  BEAM 
traces,  the  classical  spectral  parameters  arc  of  some  importance  as  weight  is  given  to  relative  signal 
power  in  the  low  frequency  bands  (parameters  1  and  3  in  Fig.  3).  Another  way  of  illustrating 
discriminant  performance  is  shown  in  Fig.  8,  where  the  separation  between  the  presumed  explosion 
(NE)  and  earthquake  (EO)  is  very  good  using  only  the  most  informative  discrimination  features, 
namely,  no.  10  and  no.  19. 

The  results  of  P  and  P-coda  parameter  discrimination  obtained  for  the  group  filter  (OGF)  and  beam 
traces  respectively  arc  shown  in  Table  3.  For  example,  for  Eastern  Kazakh  explosions  only  there 
were  no  event  classification  errors  when  using  optimum  OGF  discriminant  parameters  (Case  6, 
Table  3),  while  in  contrast  the  corresponding  BEAM  discriminant  erred  for  3  events  (Case  3,  Table 
3).  Including  the  explosions  from  the  7  other  test  sites  (Case  5,  Table  3)  gave  that  100%  of  the 
Semipalatinsk  and  75%  of  NE  from  the  other  test  sites  were  correctly  classified.  We  take  this  to 
imply  that  the  chosen  discriminant  parameters  retained  their  source-type  sensitivity  despite 
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considerable  differences  in  source  environments  and  propagation  paths.  Likewise,  the  two  S. 
Sinkiang  explosions  (test  site  no.  6;  event  nos.  17  and  18  in  Table  1)  have  similar  propagation 
paths  as  many  of  the  EO  used  in  our  analysis. 

The  EPMC  measure  is  not  linearly  related  to  the  number  of  misclassifications  and  the  same  applies 
to  those  events  which  correct  classifications  are  at  best  marginally  significant.  This  is  seen  from 
Figs.  6  and  8  and  with  further  details  in  Table  3.  The  parameters  no  10  and  19  derived  from  OGF 
traces  have  clearly  the  best  performance  without  any  misclassification  when  the  NE  populations  are 
limited  to  the  E.  Kazakh  test  site  (Case  6;  Table  3).  The  BEAM  parameters  produce  relatively 
many  misclassified  explosions,  which  with  one  exception  (event  26)  do  not  stem  from  the 
E.  Kazkh  test  site.  Event  26  remains  misclassified  even  when  the  explosion  population  is  limited  to 
the  Semipalatinsk  test  site.  Note  that  the  difference  in  the  classiOcation  errors  between  EPMC  and 
Jack-Knife  is  because  the  latter  estimates  the  miscIassiOcation  probability  for  fixed  learning  sets 
(conditional  probability),  while  EPMC  takes  account  of  randomness  in  the  learning  populations. 

We  further  remark  that  the  Jack-Knife  error  estimation  is  most  often  used  in  similar  studies. 

Among  the  earthquakes,  no.  30  in  Table  2  appears  to  be  the  most  problematic,  since  it  is 
consistently  misclassified  both  by  BEAM  and  OGF  parameters,  except  for  one  marginal  rating  in 
the  latter  case  (Case  6:  Table  3).  For  other  earthquakes  being  misclassified  or  given  a  marginal 
rating,  there  is  little  overlap  between  the  respective  OGF  and  BEAM  discriminators.  In  Fig.  9  a 
number  of  events  arc  displayed  for  which  the  discriminators  were  less  effective  or  failed.  We 
would  specifically  comment  on  these  events  in  the  next  section. 

DISCUSSION 

Firstly,  we  would  comment  briefly  on  travel  path  and  receiver-end  structural  features  which  may 
affect  the  performance  of  the  discriminators  used.  The  P  travel  time  curve  has  several  branches  in 
the  15-38  deg  distance  range  caused  by  the  so-called  400  km  and  650  km  velocity  discontinuities. 
Some  of  the  corresponding  secondary  arrivals  may  occasionally  be  very  prominent  as  demonstrated 
by  King  and  Calcagnilc  (1974)  on  the  basis  of  NORSAR  recordings  of  USSR  nuclear  explosions. 
For  the  E.  Kazakh  (Semipalatinsk)  explosions  such  secondary  phases  would  not  affect  the  coda  due 
to  a  distance  around  38  deg.  Some  of  the  explosions  are  Ored  in  basin  areas  and  the  large 
impedance  contrasts  between  sedimentary  strata  and  the  crystalline  basement  may  cause  strong 
reverberations  and  thus  affect  the  P  coda.  A  typical  example  here  is  the  Novaya  Zemlya  test  site. 
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Deslruciivc  P-pP  interference  may  affect  the  relative  spectral  power  below  1  Hz,  but  according  to 
Vergino  (1989a,b)  some  of  the  USSR  explosions  are  shallow,  that  is,  have  focal  depths  in  the 
range  150-250  m  only.  Finally,  ripple  firing  or  sequential  detonations  with  5  to  10  sec  intervals 
cannot  always  be  ruled  out.  In  the  latter  case,  coda  based  discriminators  arc  likely  to  fail. 

From  Table  3  we  have  that  the  poorest  performance  takes  place  for  the  combination  of  all  events 
and  all  discrimination  parameters  cither  derived  from  OGF  or  BEAM  traces.  Rather  surprisingly, 
the  number  of  marginal  events  seemingly  is  independent  of  event  population  and  discrimination 
parameters  used.  Regarding  the  misclassified  NE  populations,  most  of  the  events  here  are  non- 
S^;,jpaiatinsk  (E.  Kazakh),  the  exception  being  event  26  (BEAM)  and  events  22  and  5  (OGF  - 
Case  4:  Table  3).  Under  optimum  conditions  (Case  6:  Table  3),  event  5  is  labelled  marginal.  From 
the  visual  inspections  of  traces  (Fig.  9),  some  of  the  failures  are  rather  obvious,  but  this  is  not  the 
case  for  NE  26  (BEAM  -  Cases  1-3:  Table  3).  We  note  that  the  two  problematic  E.  Kazakh  events, 
no  26  and  5,  arc  assigned  ihc  smallest  m^  values  (Table  1). 

Also  some  of  the  presumed  earthquakes  arc  problematic,  in  particular  EQ  17,  22  and  32  (Cases  3 
and  5:  Table  3).  The  BEAM  and  OGF  displays  in  Fig.  9  give  the  visual  impressions  that  EO  17 
and  32  exhibit  P  waves  rather  typical  for  explosions.  A  common  feature  of  the  3  events  is  their 
northernmost  latitude  of  47°  N,  while  corresponding  longitudes  are  around  83.3°,  89.7°  and  73.6°  E 
(Table  2).  Also,  EQ  32  is  weak  (m^  -  4.3)  and  besides  look  place  in  an  area  with  no  previous 
seismic  activity.  The  ISC  bulletins  give  a  normal  focal  depth  of  33  km,  and  that  Garm  is  the  only 
USSR  reporting  station. 

Note,  most  discrimination  procedures  tend  to  favor  a  relatively  better  classification  performance  for 
earthquakes  as  also  seen  in  Table  3.  The  reason  is  that  explosion  discriminants  are  tied  to  relatively 
weak  coda  excitation  and/or  P-signal  power  deficiency  in  low  frequency  bands.  For  decreasing 
SNR  the  noise  contribution  would  be  more  relatively  important  for  explosions  and  the  net  effect  is 
that  all  weak  events  would  "look"  like  earthquakes.  In  the  latter  case  it  seems  reasonable  to  select 
for  learning  sets  weak  events  only  and  then  take  advantage  of  the  OGF-technique  for  noise 
suppression. 

From  the  above  we  have  that  the  performance  is  best  when  the  two  learning  populations  are 
drawn  from  roughly  the  .same  area.  It  is  al.so  clear  that  the  discriminators  derived  from  the  OGF- 
traccs  not  surprisingly  outperform  the  corresponding  BEAM  parameters.  This  is  not  surprising  for 
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the  simple  reason  that  the  signal-to-noise  ratio  is  signiHcantly  belter  on  the  OGF  traces  (e.g.,  sec 
EQ  32  in  Fig.  9).  Not  easily  explainable  is  the  fact  that  the  BEAM  parameters  are  weighted  in 
favor  of  the  low  frequency  part  of  the  P  signal,  slightly  in  contrast  to  that  for  the  OGF  parameters 

Our  discrimination  resulLs  compare  favorably  with  those  of  Tjostheim  (1978,  1981)  and  Sandvin 
and  Tj0sthein  (1978),  who  also  found  that  the  dominant  discrimination  power  was  embedded  in 
the  P-wave  coda.  Their  di.scriminator  parameters  were  lied  to  autoregressive  spectral  coefficients, 
which  in  turn  arc  related  to  our  prominent  parameter  no  19  (Fig.  3).  On  the  other  hand,  the 
mentioned  third  moment  discriminant  weights  most  heavily  the  presence  of  high  frequency  energy 
in  the  recorded  P  signal  and  was  reported  to  work  very  well  for  USSR  explosions  recorded  at  the 
Hagfors  (Sweden)  array  (Dahlman  and  Israelson,  1973).  This  contrasts  partly  with  our  Ondings  and 
also  with  the  results  of  Tjdstheim  (1981),  Evemden  (1977)  and  others.  The  rationale  here  is  that 
potential  high  frequency  discrimination  power  is  lost  (attenuated)  for  propagation  paths  exceeding  a 
few  hundred  kilometers.  As  also  demonstrated  here,  some  low  frequency  P-signal  discriminants  are 
important,  as  recently  reported  by  Taylor  and  Marshall  (1991)  as  well. 

Discriminators  developed  for  seismic  source  identifications  at  local/regional  distances  are  mainly 
tied  to  spectral  ratios  between  phases  like  Pn,  Pg  and  Lg  (e.g.,  see  Dysart  and  Pulli,  1990;  and 
Dowla  et  al,  1990;  sec  also  Su  et  al,  1991).  In  the  latter  study,  the  Pg  phase  was  associated  with  a 
time  window  tied  to  the  group  velocities  of  6.0  and  5.0  km/sec  -  to  us  this  is  very  close  to  a  Pn 
coda  definition.  At  these  distance  ranges,  the  Lg  parameters  have  proved  to  be  important. 

However,  a  potential  problem  with  the  Lg  pha.se  is  its  potential  blockage  for  propagation  paths 
across  sedimentary  basins  and  similar  kinds  of  crustal  heterogeneities  (Baumgardt,  1990).  The 
physical  explanation  for  the  reported  successes  in  source  identifications  at  local  and  regional 
distances  is  attributed  to  stronger  shear  wave  excitation  for  EQ-type  sources  (strong  Lg),  which 
also  would  imply  stronger  P  coda  at  any  distance. 

Finally,  we  want  to  comment  on  various  discrimination  techniques  in  use.  At  the  outset  of  this 
study  the  "neural  network"  approach  was  not  contemplated,  and  besides  Dysart  and  Pulli  (1990) 
and  Dowla  et  al  (1990)  reported  that  the  Fisher  discriminant  had  equivalent  performances.  The 
reason  for  this  is  that  the  EQ  and  NE  populations  are  essentially  linearly  separated.  It  is  here  meant 
that  the  feature  parameters  derived  from  the  respective  event  populations  are  separated  by  a  straight 
line  in  the  two  parameter  case  (Fig.  7).  We  would  here  add  that  our  approach  is  relatively  robust 
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for  small -sized  cvcni  populations  and  furihcrmorc  that  the  probability  of  misclassifiaition 
introduced  here  is  helpful  in  judging  discrimination  performances. 

Both  the  "neural  network"  and  the  Fisher  discriminators  arc  constructed  for  multidimensional 
classification  problems,  in  our  case  u  priori  known  NE  and  EO  event  populations.  Strictly 
speaking,  this  is  not  ihe  ease  lor  the  non-Scmipalalinsk  arca.s,  so  not  unexpectedly  the 
discrimination  pcriormanccs  were  poorer  for  the  other  USSR  tost  site.s.  An  alternative  strategy  for 
such  cases,  which  could  he  highly  relevant  in  many  contexts,  is  to  consider  singular  learning 
samples  of  cither  explosions  or  earthquakes  (e.g.,  sec  Tsvang  ct  al,  1986).  The  problem  is  related 
to  the  task  of  finding  the  most  powerful  di.scrimination  parameters  which  clearly  would  affect  to 
some  extent  the  source-receiver  wavepath  and  the  local  structure  in  the  source  area.  Our  approach, 
tied  to  testing  an  initial  large  cla.ss  of  discriminators  and  then  estimating  which  ones  are  most 
powerful,  appears  to  be  highly  efficient  in  practice. 

CONCLUSION 

In  this  study  we  have  presented  a  comprehensive  seismic  source  discrimination  .scheme  for 
Icleseismic  events,  rite  observational  data  were  ba.scd  on  USSR  events  as  recorded  by  the 
NORLSS  (Norway)  array 

Th''  learning  population  comprised  44  nuclear  explosions  (NE)  and  35  earthquakes  (EO)  from  the 
general  E.  Kazakh  area.  An  explosion  subset  comprksed  32  events  from  Scmipalatinsk.  The 
di.scriminant  features  used  in  analysis  were  extracted  both  from  the  conventional  array  beam  traces 
and  from  lx;am  traces  derived  by  ()(iF,  where  low-frequency  (below  2  Hz)  coherent  noi.se  is 
suppre,ssed.  Major  results  were  as  follows; 

I  A  total  ol  19  potential  discrimination  features  wcie  consideied.  These  were  mainly  lied  to 

normalized  spectral  power  m  8  frequency  bands  for  both  P-  and  P-coda  waves  of  durations 
6  sec  and  30  .sec,  respectively  Also,  the  clas.sical  P-complexity  discriminant  was 
incorporated. 

2.  The  most  cllectixe  ilisciiminaiion  leatures  weie  the  complexity  one,  and  relatively  low  P- 

signal  Ifcqucncies.  The  best  cla.ssiliatlion  perlbrmancc  was  rrirtained  by  using  2  to  4 
features  which  in  turn  were  .selected  on  the  basis  of  the  extracted  measure  of  maximum 
misclassificalion  probability. 


V  The  rcntiitus  extracted  Iroin  the  0(J1  event  traces  had  a  sliglitly  l)etter  [X'rlorinance  than 

those  derived  from  the  conventional  event  beam  traces.  I’lie  corresponding  optimum  Icature 
subsets  were  somewhat  dilTerent. 

4.  When  the  explosion  population  was  restricted  to  the  Semipalaiinsk  area,  a  complete  event 
classification  was  obtained  using  OGl'-derived  features.  However,  3  events  here  were  rated 
of  marginal  significanee. 

.S.  Inclusion  of  all  the  Nli  Irom  the  other  7  test  sites  resulted  in  a  slight  increase  of  6.5  [kt 
cent  in  the  classification  error.  This  implies  that  the  chosen  discriminant  parameters  retain 
their  source-type  sensitivity  despite  considerable  differences  in  source  environments  and 
propagation  paths.  Iakcwi.se,  our  dllscriminants  were  able  to  separate  S.  Sinkiang 
explosions  from  nearby  earthquakes  having  .similar  propagation  paths. 

6.  A  remaining  problem  is  that  of  designing  discriminants  for  areas  where  available  events  are 
essentially  limited  to  cither  earthquake  or  explosion  populations. 
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lAHIX  1 


I'tcsuii'tJ  lixplosioiis 


No 

Dale 

Rcpioii 

I  jti. 

IjOii 

Mb 

Mi 

1 

1934 

12 

28 

East.Ka:.  (3) 

49.84 

78.75 

4.0 

4.1 

2 

1985 

02 

10 

East.Xaz.  (5) 

49.87 

78.31 

5.9 

4.4 

3 

1935 

04 

25 

East.Ea:.  (5; 

49.92 

78.97 

5.9 

5.2 

1985 

04 

30 

Eist.Ku:.  (5) 

49.64 

73.69 

4.0 

4.2 

5 

1985 

07 

11 

East.Kaz.  (5) 

49.90 

78.80 

3.5 

— 

6 

1985 

07 

18 

Vest. Russia  (7) 

45.97 

40.64 

5.0 

3.7 

7 

1935 

07 

25 

East.Kaz.  (5) 

49.89 

78.15 

5.0 

4.0 

3 

1988 

09 

14 

East.Kaz.  (5) 

49.82 

78.79 

4.1 

4.4 

9 

1988 

11 

12 

East.Kaz.  (3) 

50.05 

78.99 

5.2 

— 

10 

1983 

11 

23 

East.Kaz.  (3) 

49.78 

78.14 

5.3 

— 

11 

1987 

02 

24 

East.Kaz.  (5) 

49.81 

78.14 

5.4 

— 

12 

1987 

03 

12 

East.Kaz.  (3) 

49.93 

78.78 

5.4 

3.9 

13 

1987 

04 

03 

East.Kaz.  (5) 

49.90 

78.80 

6.2 

4.7 

If 

1987 

04 

17 

East.Kaz.  (3) 

49.83 

78.69 

6.0 

4,3 

13 

1987 

04 

19 

Ural. Haunt.  (8) 

40.78 

54.22 

4.3 

— 

16 

1987 

03 

04 

East.Kaz.  (5) 

49.77 

-78.09 

5.5 

— 

17 

1937 

04 

03 

S.SinIrianq  (4) 

41.58 

88.73 

6.3 

18 

1987 

04 

04 

East.Kaz.  (3) 

49.84 

78.14 

5.4 

— 

19 

1987 

07 

04 

(Cent, Siberia  (3) 

41.49 

112.78 

■5.2 

‘  — 

20 

1987 

07 

17 

East.Kaz.  (3| 

49.78 

78.12 

5.8 

1*^ 

21 

1987 

07 

24 

;Cent.Siberia  (3) 

61.46 

112.72 

5.1 

/  V 

22 

1987 

08  02- 

'Eastikaz.  (3) 

,49.84 

:  78i88 

5.9 

3.8 

23 

1987 

08 

.02 

N.lMlya  (1) 

73.3t‘ 

54.71 

5.8 

5.4 

24 

1987 

08 

12 

Cent. Siberia  (3) 

61.42 

112.71 

5.0 

— 

23 

1987 

10 

03 

■  Hest.taz.  (4) 

47.63 

36.21 

5.2 

»  x  -  ■ 

24 

1987 

10 

14 

East.Kaz.  (3) 

49.78 

78.24 

4.4 

27 

1987 

11 

15 

East.Kaz.  15) 

49.87 

78.79 

4.0 

4.8 

28 

1987 

12 

13 

East.Kaz.  (3) 

49.95 

78.85 

4.1 

4.5 

29 

1987 

12 

20 

East.Kaz.  (3) 

49.75 

78.02 

4.8 

— 

30 

1987 

12 

27 

East.Kaz.  (SI 

49.83 

78.74 

4.1 

4.5 

31 

1983 

02 

04 

East.Kaz.  (3) 

49.84 

77.96 

4.8 

— 

32 

1988 

02 

13 

East.Kaz.  (3) 

49.92 

78.90 

6.0 

4.5 

33 

1988 

04 

03 

*East'.Kazl  (3) 

'49^88 

78.96 

5.9 

— 

34 

1988 

03 

04 

East.Kaz.  (3) 

49.91 

78.72 

4.2 

— 

35 

1988 

05 

07 

N.ZesIya  (1) 

73.35 

54.24 

5.5 

3.8 

34 

1988. 

04 

14 

East.Kaz.  (31 

50.02 

78.98 

4.9 

4.1 

37 

1988 

08 

22 

Vest. Siberia  (2) 

46.28 

78.55 

5.3 

— 

38 

1988 

09 

29 

S.Sintianq  (4| 

41.82 

88.25 

4.8 

— 

39 

1983 

10 

18 

East.Kaz.  (3) 

49.84 

78.12 

4.9 

40 

1998 

12 

04 

N.Zetlya  (1) 

75.36 

55.07 

5.9 

— 

41 

1988 

12 

17 

East.Kaz.  (31 

49.88 

78.92 

5.9 

4.5 

42 

1989 

02 

12 

East.Kaz.  (3| 

49.89 

78.75 

5.9 

4.4 

43 

1989 

02 

17 

East.Kaz.  (5) 

49.87 

78.17 

5.0 

44 

1989 

22 

05 

East.Kaz.  (51 

49.91 

78.86 

6.0 

4.5 

Table  I.  Listing  of  Ibe  presumeti  Mndergsound  nuclear  explosions  (NE)  recorded  at  NORESS  and 
used  in  our  analysis.  The  region  numbeis  (in  brackets)  refer  to  the  spedfic  USSR  lest  sites  shown 
in  Fig  1  The  focal  parantetets  ire  taken  from  the  PDE  and  iSC  bulletins.  Note,  Event  11  (Ural 
Mountains)  ii  the  lecond  explosloo  on  19  April  with  origin  lime  at  04.04.SS.6.  The  list  is 
incooiplele  as  NORESS  data  limn  5  preaumod  NEs  for  uaifous  techalcal  reasons  were  not  available 
10  us.  TTie  events  are:  15  Jus  1985  (48E9N.  78E8EX  7  Nov  1985  (49.95N,  7&83E),  19  Apr  1987 
(60.67N,  56J0E).  20  Jon  1987  («9.9Wi,  78.73E)  and  6  Sep  1988  (61.22N.  48.04E).  Another  ‘tost' 
event,  possibly  an  NE,  Is  that  of  18  Sep  1987  (50.18N,  78.02E). 
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1  ABI.l-  I 


Presumed  Earthquakes. 


No 

Dale 

Lat. 

Lon. 

Depth 

Mb 

Ms 

1 

1984  11  02 

42.00 

84.06 

33 

4,5 

•> 

4. 

1985  02  03 

42.06 

84.34 

- 

4.5 

— 

3 

1985  03  24 

42.06 

77.62 

3 

4.3 

— 

4 

1985  04  16 

42.26 

82.24 

33(3) 

4.5 

— 

5 

1985  06  02 

43.80 

85.65 

21(45) 

4.9 

4.2 

6 

1985  07  16 

42.22 

82.36 

45(22) 

4.9 

^2 

7 

1985  08  14 

42.13 

82.44 

33 

4.1 

— 

8 

1985  08  23 

42.65 

74.50 

3 

4.3 

— 

9 

1986  05  30 

43.23 

87.82 

33 

4.6 

— 

10 

1986.06  12 

43.67 

87.33 

33 

4.8 

. — . 

11 

1986.07  03 

43.91 

84.69 

.  33 

4.4. 

12 

4986^07.'17 

43.32 

77.85 

33 

-4.5 

13 

;i9^.07  21 

44.66 

79.50 

33 

■4  .'6 

H 

1986*07  24 

43.79 

87.25 

33 

4.5 

- 

IS 

1986  10  04 

42.39 

84.63 

‘33 

4.0  . 

' - : 

16 

,1986  11^0 

42.04 

84.39 

50 

4.6 

- ■ 

17 

Xm'Xl  14 

47.31 

83.31 

33 

5.0 

18 

'iw’ooM 

42,45 

79.95 

33 

19 

il987^5'l0 

44.28 

79.74 

33 

4.5 

20 

1987  05  26 

42.92 

78.06 

20 

4.6 

21 

1987  08,22 

43.81 

85.29 

58 

4.4 

— . 

22 

1987  09' 18 

47.01 

89.65 

33 

5.3 

4.8 

23 

1987  09  20 

42.91 

77.61 

41 

4.6 

4.2 

24 

1987  10  06 

43.43 

88.54 

32 

4.8 

4.0 

25 

1987  10  16 

44.20 

82.84 

56 

4.7 

— 

26 

•1988  02^08 

43.73 

83.76 

10 

.4.3 

- - T. 

27 

Xm  03*13 

42.18 

75.44 

33 

4.5 

- - 

28 

•1988  03;15 

42.21 

75.50 

33 

4.5 

29 

1988  03 '25 

44.70 

79.60 

33 

4.5 

— 

30 

1988  05  25 

42.01 

85.69 

22 

5.2 

— 

31 

1988  06  17 

42.97 

77.50 

24 

5.3 

5.3 

32 

1988  09  27 

46.80 

73.59 

5 

4.3 

— 

33 

1988  11  15 

42.01 

89.29 

33 

5.0 

4.3 

34 

1989  03  05 

42.51 

74.63 

33 

5.3 

4.1 

35 

1989  05  08 

44.83 

79.92 

33 

4.7 

— 

Table  2.  Listing  of  presumed  earthquakes  recorded  at  NORESS  and  used  in  analysis.  Focal 
parametets  are  taken  from  (SC  and  POE.  The  numbers  in  brackets  correspond  to  PDE  depths  when 
different  from  ISC  The  choice  of  this  EQ  region  (insert  in  Fig.  1)  reflects  the  fret  that  only  2  test 
sites  in  Eurasia,  that  is,  Semipalatinsk  and  S.  Sinkiang,  arc  reasonably  close  to  areas  of  high 
seismicity.  The  list  is  inoomplele  as  NORESS  recordings  from  14  presumed  earthquakes  were  not 
available  to  us  for  various  technical  reasons. 
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TABLB  3. 


DISCRIMINATION  RESUL'I'S 


Case  Data  Populations  Class,  errors  Qass.  Misclass.  Marginal 

Rcclass.  J. knife  features  events  events 


1 

BEAM  All  events 

7 

9 

All 

NE:6,25,26,35, 

3738; 

NE:15; 

EO:5,13,18; 

EQ:22,2335; 

2 

BEAM  All  events 

3 

7 

1,193.13 

NE:6,23,26,35, 

37,38; 

NE:5,40; 

EQ:23; 

EQ:13,31; 

3 

BEAM  All  EQ 

1 

3 

19,13 

NE:26; 

Kazakh  NE 

EQ:22,32; 

EQ:23; 

4 

OGF  All  events 

5 

8 

All 

NE:5,6,22,23 

3738; 

NE:25; 

EQ:17,32; 

EQ: 13,22,23; 

5 

OGF  All  events 

3 

5 

19,10,15, 

NE:6,26,28 

NE:23; 

13 

EQ:17,32 

6 

OGr  All  EQ 

0 

0 

19,10 

— 

NE:5,29,34; 

Kazakh  NE 

EQ:  13,32; 

Table  3.  Discrimination  results  for  the  Eurasian  events  listed  in  Tables  1  and  2,  and  recorded  at 
the  NORESS  array.  Qassification  features  arc  extracted  from  ordinary  beam  (BEAM)  traces  and 
optimum  group  Gitered  (OGF)  traces,  where  coherent  noise  has  been  suppressed.  The  explosion 
(NE)  population  was  divided  in  two  parts:  all  explosions  or  only  those  at  the  E.  Kazakh  test  site 
(see  Fig.  1).  QassiGcation  (class.)  enots  are  given  for  two  cases:  ReclassiGcation  where  the  event 
tested  was  part  of  the  learning  population  and  Jack-knife,  where  test  event  and  learning  population 
were  independent.  The  "Reclass.”  results  are  always  somewhat  optimistic.  ClassiGcation  features 
are  detailed  in  Fig.  3.  Misclassified  events  and  Marginal  events  refer  to  the  numbering  in  Table  1 
and  2.  The  most  problematic  event  is  EQ  32,  that  is,  a  presumed  earthquake  of  27  Sep  1988 
(46.80N,  783.59E),  which  also  has  a  visual  appearance  of  being  an  explosion.  Waveforms  for  a 
representative  number  of  events  being  either  misclassified  or  rated  marginal  are  shown  in  Fig.  9. 
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classification  study.  The  test  sites  arc  N.  Zcmlya  (1);  W.  Siberia  (2);  Cent.  Siberia  (3);  W.  Kazakh  (4);  E. 
Kazakh  (5);  S.  Sinkiang  (6);  W.  Russia  (7);  and  Ural  Mountains  (8).  Focal  parameter  details  in  Tables  1  and  2. 


POWER  SPECTRA  (8  bands) 


o>  *o  O 
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here  is  aa  earthquake  -  no.  21  in  Table  2. 


PARAMETERS 
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bands  used  are  shown  to  the  left.  Note,  all  spectral  band  values  were  normalized  for  each  event  with 
respect  to  Amax.  In  the  calculation,  the  logarithmic  feature  values  were  used. 


P-spectra 


44  explosions 


Fig.  4.  NomMitze  power  spectra  (signal-noise)  for  all  events  osed  in  analysis  (Table  1  and  2).  'Peaks' 
correspond  to  the  center  frequency  of  the  eight  spectral  bands  shown  in  Fig.  2.  Seemingly,  the  largest 
spectral  differences  between  the  explosion  and  earthquake  populations  are  in  the  range  O.S  -  l.S  Hz. 


pakamcicks  SEi.rcTioN  (ocr) 


Ranked  parameters  (o-all  events;  +-only  Semipol.) 


PARAMETERS  SELECTION  (BEAM) 


Ranked  parameters  (o-all  events;  -t — only  Semipol.) 

Eij!-  6.  The  relative  performance  of  feature  parameters  and  combinations  hereof  as  a  function  of  the 
Expected  Probability  of  Misclassiflcation  (EPMC)  as  denned  in  the  text.  OOF  refers  to  optimuni  group 
filtered  traces  (removal  of  coherent  noise),  while  BEAM  refers  to  ordinary  beam  traces.  Also,  the  upper 
curves  refer  to  all  events,  while  the  lower  ones  refer  to  the  subset  of  the  Semipalatinsk  or  E.  Kazakh  test 
site  explosions  (Fig.  1).  The  minimum  in  the  EPMC  parameter  implies  that  a  combination  of  2  or  4 
classiftcation  features  is  optimum,  as  inclusion  of  additional  feature  parameters  actually  decreases  the 
performance  as  discussed  in  the  text. 


U9 


Log  OGF  (optimal  parameters:  IO-t-19) 


Fig.  7.  Diagram  for  Uie  feature  parameters  10  and  19  as  calculated  from  CXJF  traces  for  all  earthquakes 
(Table  2),  but  the  explosion  population  was  limited  to  the  E.  Kazakh  test  site  (Table  1).  The  two  event 
populations  are  well  separated,  but  a  better  display  is  given  in  Fig.  8b,  while  Case  6  in  Table  3  actually 
details  marginal  events. 
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HI  AM  DISCRIMINATION 


.^bl’0  +  44NC,all  parctiietcr  y 


35EQ  +  44NE,  opiimal  parameiers  (1,19,3,13) 


Only  E.Kazach  (32NE+35EQ).  optimal  parameters  (19,13) 


Evcnis  numberj  ranked  due  fo  UDF  values 


A-iTvsciassifted  euplosions  e-miscJassiried  earthquakes 

Fig.  8a.  Event  discrimination  with  feature  parameters  extracted  from  the  conventional  beam  traces.  The 
l.ni-  (linear  diseriniination  fnm  uon)  is  defined  in  the  text  and  Jack-Knife  refers  to  the  specific  way  of 
testing  individual  events  Mid<'|K’ndciitly  against  the  learning  jxipulation  The  cvenLs  arc  ranked  relative  to 
their  EDI-'  values,  and  thu.s  not  the  way  done  in  Tables  1  and  2.  Events  with  LDF  values  within  ±  0  2 
units  (stippled  lines)  arc  somewhat  arbitrarily  defined  as  marginal  events.  The  3  combinations  of  event 
populations  and  optimum  p.iiamcier  combinations  coincide  with  Case  1*3  in  Table  3. 
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OGF  DISCRIMINATION 


35EQ  +  44NC,all  porarneters 


3 

^  2 

« 

C  1 
I 


o 

-1 

L.- 

o 

-  -2 
-3 


35EQ+44NE,  optimal  parameters  (19,10,15,13) 


Only  E.Kozach  (32NE+35EQ),  optimal  parameters  (19,10) 


A-m>sciass<fi«d  explosions 


•-misdassified  earthquakes 


Fig.  8b.  Event  discrimination  with  feature  parametcis  extracted  from  the  OGF  traces  (optimum  group 
filtered).  Caption  otherwise  as  tor  Fig,  8a. 
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